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Primordial gravitational waves provide a very important stochastic background that could be 
detected soon with interferometric gravitational wave antennas or indirectly via the induced patterns 
in the polarization anisotropics of the cosmic microwave background. The detection of these waves 
will open a new window into the early Universe, and therefore it is important to characterize in 
detail all possible sources of primordial gravitational waves. In this paper we develop theoretical 
and numerical methods to study the production of gravitational waves from out-of-equilibrium 
gauge fields at preheating. We then consider models of preheating after hybrid inflation, where the 
symmetry breaking field is charged under a local U(l) symmetry. We analyze in detail the dynamics 
of the system in both momentum and configuration space. We show that gauge fields leave specific 
imprints in the resulting gravitational wave spectra, mainly through the appearence of new peaks at 
characteristic frequencies that are related to the mass scales in the problem. We also show how these 
new features in the spectra correlate with string-like spatial configurations in both the Higgs and 
gauge fields that arise due to the appearance of topological winding numbers of the Higgs around 
Nielsen-Olesen strings. We study in detail the time evolution of the spectrum of gauge fields and 
gravitational waves as these strings evolve and decay before entering a turbulent regime where the 
gravitational wave energy density saturates. 
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I. INTRODUCTION 

Gravitational waves (GW) are a robust prediction of General Relativity PQ. They correspond to ripples in space- 
time that travel at the speed of light, and are typically produced whenever an astronomically large body of mass moves 
at relativistic speeds like in astrophysical binary systems, or whenever large density contrast waves collide against 
each other, like in early universe phase transitions. The change in the orbital period of a binary pulsar known as 
PSR 1913+16 was used by Hulse and Taylor [2] to obtain indirect evidence of their existence. Although gravitational 
radiation has not been directly detected yet, it is expected that the present universe should be permeated by a diffuse 
background of GW of either astrophysical or cosmological origin [1 . Astrophysical sources, like the gravitational 
collapse of supernovae or the neutron star and black hole binaries' coalescence, also produce a stochastic gravitational 
wave background (GWB) which comes from unresolved point sources. On the other hand, among the backgrounds 
of cosmological origin, we find the approximately scale-invariant background produced during inflation [5], or the 
GWB generated at hypothetical early universe thermal phase transitions [3HI] , from relativistic motions of turbulent 
plasmas [5] or from the decay of cosmic strings [1 . Fortunately, these backgrounds have very different spectral 
shapes and amplitudes that might, in the future, allow gravitational wave observatories like the Laser Interferometer 
Gravitational Wave Observatory (LIGO) 9 , the Laser Interferometer Space Antenna (LISA) [TO], the Big Bang 
Observer (BBO) QT] or the Decihertz Interferometer Gravitational Wave Observatory (DECIGO) [IS], to disentangle 
their origin. Unfortunately, due to the weakness of gravity, this task will be extremely difficult, requiring a very 
high accuracy in order to distinguish one background from another. It is thus important to characterize the many 
different sources of GW as best as possible. Recently there has been a significant improvement in the sensitivity 
of laser interferometers to a cosmological background of GW and there are now limits on the amplitude of this 
background that are just below the BBN [13] and CMB bounds [H] as recently reported by the joint LIGO and 
VIRGO collaborations [H]. 

One cosmological GW background that is very well motivated by other observations is the approximately scale- 
invariant spectrum of GW produced from quantum fluctuations during inflation [3], This spectrum extends over a 
very wide frequency range and its amplitude is directly related to the energy density during inflation. These GW may 
be detected indirectly by forthcoming CMB experiments through their effect on the B-mode polarization of CMB 
anisotropies [16] . or even directly in the longer term by interferometric experiments, if inflation occurs at sufficiently 
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high energy scales. However, if inflation occurs at lower energies, as is the case in many models motivated by high- 
energy physics, GW from inflation would have an amplitude that is too low to be observable. On the other hand, 
in case that a global phase transition took place after the end of inflation, then a scale-invariant GW background is 
also generated [TTrtl9| . The detailed calculations show that the amplitude of this new GW background is two orders 
of magnitude greater than that expected from inflation for the same energy scale |19j . and thus might be detected 
directly by the future GW detectors or in the B-mode polarization of the CMB [20]. Finally, another source of GW 
that may be relevant for interferometric experiments and whose study will be our main target here, is provided by the 
violent period following the end of inflation. Indeed, in many models, the inflaton decays in an explosive and highly 
inhomogeneous way, in the process of preheating [21j . The particular mechanism responsible for preheating is model- 
dependent, but it is generally dominated by a non-perturbative production of Bose fields with very high occupations 
numbers, far from thermal equilibrium. This leads to a second, longer stage characterized by turbulent-like interactions 
between classical waves, before the system eventually reaches a thermal state [55]. The large, time-dependent fields' 
inhomogeneities produced by preheating source a stochastic GW background [2"3rt3T] . When redshifted until today, 
this background may fall in the frequency range accessible by interferometric experiments if inflation and preheating 
occur at sufficiently low energy scales, providing an alternative to test inflation with GW. In addition, GW from 
preheating carry crucial information about the mostly unkown post-inflationary dynamics and, because the details 
of preheating depend very much on the model, they could be used in the future to discriminate between different 
inflationary models. 

Gravitational waves from preheating have been intensively studied by different groups |23H31| . Different numerical 
methods have been used in early works, see [3DJ for a critical comparison, but the later results obtained by different 
groups agree well with each others [26] |29TI31| . Two main classes of models have been studied so far: preheating 
after chaotic inflation and preheating after hybrid inflation. In the first case, preheating proceeds via parametric 
resonance [21] • Because inflation and preheating occur at high energy scales in these models, the resulting GW have 
a typical frequency today in the range 10 6 - 10 9 Hz, which is too high for the signal to be observable by currently 
available experiments. On the other hand, preheating after hybrid inflation may occur at much lower energy scales and 
the resulting GW may fall in a frequency range below 10 3 Hz, which is accessible by high-sensitivity interferometric 
experiments. GW from preheating after hybrid inflation [35J [33J were first studied in Ref. [23], in the framework of 
parametric resonance studied in |23) . It was later understood that hybrid inflation models preheat in an even more 
violent way, due to the tachyonic amplification of fluctuations of the symmetry breaking field when the fields roll 
towards the true minumum of the potential, a process called tachyonic preheating 34J. Gravitational waves from 
tachyonic preheating after hybrid inflation were first studied in Ref. [2"5H5D] and then fully explored in [3T] , where the 
regions of the parameter space that may lead to an observable signal were determined. Finally, it is also worth noting 
that the methods developed to study GW production from preheating may be applied to other out-of-equilibrium 
sources in the early universe. One example that is particularly interesting from the perspective of ground based 
interferometers is the GW background [35J produced by the non-perturbative decay of flat direction condensates in 
supersymmetric theories [36] . Another example is the evolution of unstable domain walls [39] . 

One common point in these works has been to focus on models involving only scalar fields 2 . However, in realistic 
models, gauge (vector) fields may also be copiously produced during preheating [32 [3U [35] 3 . These lead to new 
terms in the anisotropic stress sourcing GW, in addition to the gradients of the scalar fields, and they may play an 
important role in GW production. Indeed, the numerical simulations [25, 29-31 indicate that scalar fields do not lead 
to a significant production of GW during the turbulent evolution towards thermal equilibrium after preheating. This 
was demonstrated in [3D] where it was shown that massless gauge fields may change this result. Out-of-equilibrium 
gauge fields are of course also ubiquitous in other sources of GW, such as thermal phase transitions and local topological 
defects [1] . Moreover, tachyonic preheating [331 H31 S3] could be responsible for copious production of dark matter 
particles [35], lepto and baryogenesis [33JI3BH3HI, topological defects [33] |33J, primordial magnetic fields [22J, etc, 
whose observational consequences could help put more stringent constraints on the period immediately after inflation 
responsible for the reheating of the Universe. 

The main purpose of this paper is to study the GW background produced by physical models with scalar and 
gauge fields at preheating. We will develop numerical methods and derive theoretical results which can be applied 
to out-of-equilibrium gauge fields in general. We will then consider models of preheating after hybrid inflation where 
the symmetry breaking field is complex and coupled to a £7(1) gauge field. As we will see, a major consequence of 



1 GW may also be produced from Q-ball fragmentation [37| . see however |38] . 

2 GW production from a decaying tachyon field coupled to a gauge field has been studied numerically in Ref. 40 , although lattice gauge 
techniques do not seem to be used in this work and the transverse-traceless part of GW does not seem to be correctly extracted. 

3 Similarly, flat direction condensates in the MSSM are charged under the gauge symmetries and gauge fields are produced by their 
non-perturbative decay |36| . 
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the gauge field is to introduce new characteristic scales which are inherited by the GW spectra. This prompted us 
to study in detail the dynamics of the system in both Fourier and position spaces. We found in particular that a 
crucial role is played by the dynamics of cosmic string configurations of the gauge and scalar fields. In order to probe 
these different scales in the simulations, we had also to develop a lattice calculation of GW with gauge fields that is 
accurate up to second order in the lattice spacing and timestep. 

The paper is divided as follows. In section [TTJ we specify the model of hybrid preheating embeded in a gauge 
framework that we will study in this paper, briefly reviewing the fields' dynamics after inflation. Section |III| is 
dedicated to theoretical perspectives on the effects of gauge fields on GW production from preheating. In section |IV] 
we present our numerical method to compute the GW spectra produced by scalar and gauge fields on the lattice. In 
section [Vj we study in detail the dynamics of an abelian-higgs model of tachyonic preheating in Fourier space and 
the resulting GW spectra. Section |VI| is dedicated to the sudy of spatial configurations in position space of both the 
sources and the GW themselves. We conclude in section [VTT| with a summary of our results and directions for future 
works. In appendix [Aj we derive a no-go result for abelian models discussed in the main text. The details of our 
lattice calculation are given in appendix |B} 

II. ABELIAN-HIGGS PREHEATING AFTER HYBRID INFLATION 

Hybrid inflation [32j is a class of inflationary models based on particle physics and in particular on spontaneous 
symmetry breaking by a scalar field. This field is coupled to a singlet scalar, the inflaton, and triggers the end of 
inflation when its mass-squared goes from positive to negative due to its coupling to the rolling inflaton. Let us 
consider a generic hybrid inflation model, as described by the potential 

v{v, X ) = - (M 2 - v 2 ) 2 + i 5 y M 2 + y infl ( x ) , (i) 

where g 2 is the strenght of the coupling of the inflaton \ to a symmetry breaking field ip, with self-coupling A and 
vacuum expectation value (VEV) v. A relevant feature of these models is that a small A is not a prerequisite in order 
to generate the observed CMB anisotropics, while the scale of inflation can be choosen to range from GUT scales 
(~ 10 16 GeV) all the way down to GeV scales. Depending on the particular model considered, the purely inflaton 
part of the potential 14 n fi(x) can take different forms, see e.g. Ref. [SO], and an apropiate choice of the parameters 
makes these models completely compatible with CMB constraints. 4 

We want to embed this setup into a gauge- invariant framework, such that the symmetry breaking field <p is coupled to 
the corresponding gauge fields. In general, the group of gauge symmetry could be abelian or non-abelian (or a product 
of both) and ip could even be the Higgs field of the Standard Model with gauge symmetry SU(3) x SU(2) x U(l). In the 
very early universe the gauge group could for instance contain multiple U(l)'s, as described in some compactifications 
of string theory. In the numerical simulations, however, we will restrict ourselves to a complex field ip = (ipi + iip2) 
coupled to a single gauge field An, with a U(l) symmetry. Nevertheless we will also consider other gauge groups in 
section |Ul| where we derive theoretical expectations about the production of GW from gauge fields. 

For convenience, independently of the physical origin of ip, we will often refer to this field simply as the "Higgs" 
of the model. In the case of a U(l) symmetry, we will consider the Abelian-Higgs model coupled to an inflaton, as 
described by 

S = -J d^x^^F^ + (D^)\D»p)+ 1 -(d, X ) 2 + V(tp,x)}, (2) 

with g the background metric determinant, — d^A u — d v A^ is the (anti-symmetric) field strength of the gauge 
field A^, and Z)„ = <9 M — ieA^ the usual covariant derivative with e the gauge coupling. After inflation, the metric is 
that of a (flat) Friedman-Robertson- Walker (FRW) background, with g^ v — diag(— 1, a 2 (i), a 2 (i), a 2 (t)) and a(t) the 
scale factor. The (classical) equations of motion can then be derived varying the action ^ as 

-d^ X + 3H X + V x = (3) 
-D^p + 2,Hp + V, r = (4) 
E i + HE i -e ijk d j B k -2elm[(D i cp)*ip] = 0, (5) 



4 For instance, the original model [32], with V infl (x) = M 2 X 2 and fi = 1.4 X 10 14 GeV, v = 3.6 X 10 16 GeV, A = 0.17 and g = 0.001, 
produces a spectrum with scalar tilt n B = 0.98 and tensor-to-scalar ratio r = 0.1, in perfect agreement at 95% c.l. with the WMAP-7yr 
data [51] . Other models based on logarithmic loop corrections to the flatness of the potential also give a negative tilt during inflation 
compatible with WMAP data. 
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where a dot denotes the time derivative, H = a/a is the Hubble rate, E % — F at and Bi — \tij k F^ k are the gauge 
invariant electric and magnetic fields and e^-fc is the completely anti-symmetric 3-tensor. Here and in what follows, 
latin indices k, ... run over the three spatial coordinates and repeated indices are summed unless stated otherwise. 
The equation associated to the time component of the gauge field, 

d k E k = -2elm[{D tp)*tp\, (6) 

is not a dynamical equation, but rather a constraint equation equivalent to the Gauss Law of classical electromag- 
netism, \7E = p, with p — — 2elm [{D^ip)* tp]. 

In this model the characteristic time scale evolution of the fields after inflation is set by the inverse of to = \/\v. It 
follows then that the expansion of the universe is negligible during preheating, since H ~ (v/M p )m -C to. For studying 
preheating it is thus sufficient to consider a flat background rj^ v — diag(— 1, 1, 1, 1). Moreover, for convenience we will 
evolve the system in the temporal gauge, i.e. fixing A — 0. The Higgs and gauge fields equations then look as 

tp - DiDiip + V )ip . = (7) 
Ai-djdjAi+didjAj = 2elm[tp* Dup] (8) 
d t A, = 2elm[tp*tp] . (9) 

The condition Aq — does not fix the gauge completely as we can still perform a gauge transformation tp — ¥ e ta tp, 
A- L — > Ai + -did, with a time-independent function a = a(x). Numerically, we will evolve the system in the temporal 
gauge and use the remaining gauge degree of freedom to set Ai = initially. The initial conditions for the electric 
field are then read from satifying initially the Gauss constraint ([9]), see below. 

The abelian Higgs model is known to give rise, upon symmetry breaking, to cosmic strings connecting Nielsen- 
Olesen vortices, where the Higgs winds n-times. The type of string (I and II) depends on the ratio of Higgs to gauge 
field mass, j3 = 2A/e 2 = to^/to^. Magnetic flux lines repel each other, while the scalar field produces an attractive 
force, and their range is controlled by the Compton wavelength of the mediating boson. Depending on which one is 
more massive we may have type I or type II superconductors. For /3 < 1 (A < e 2 /2) the penetration depth of the 
magnetic field inside the string is smaller than the coherence length of the string, and local vortices are stable for 
arbitrary winding n, as in type I superconductors, while for /3 > 1 (A > e 2 /2) we have the opposite situation, and the 
magnetic flux lines can live inside the string (vortices with n > 1 are unstable in that case). Some of the properties 
that we will encounter in the time evolution of the abelian Higgs model after symmetry breaking in fact depends on 
whether the cosmic strings that are formed are type I or II. In both cases we find cosmic strings with windings of the 
Higgs around them, and their energy density will be a strong source for gravitational wave production. We leave for 
Sections V and VI the detailed description of the rather complicated dynamics. 



A. Hybrid Preheating 

The general qualitative dynamics of the system is as follows. For \ > Xc where Xc = m /g is called the critical point, 
the Higgs field have positive mass squared and the potential has a valley at tp — 0. During inflation the expansion 
is driven by the false vacuum energy, Vq ~ Aw 4 /4, while the inflaton decreases slowly along the valley due to the 
uplifting term Vm/(x) m Q- Inflation ends either when x reaches the critical point or when the slow-roll conditions 
are violated, whichever occurs first. In both cases, when \ < Xo the effective square Higgs mass becomes negative 
and this triggers the symmetry breaking process from (\tp\ 2 ) — to the true vacuum (\tp\ 2 ) = v 2 . The time at which 
the inflaton reaches \c is called the critical time t c . Around that time, the inflaton's time evolution can be linearly 
approximated by x{t) — XcO- ~ V c m{t — t c )), where and V c is the (dimensionless) inflaton's velocity at t c , 



Xc 

Through the Higgs-inflaton coupling, the evolution of x(t) around \c induces a time dependence in the Higgs effective 
mass, m 2 = — 2V c m 3 (t — t c ), which changes from positive at t < t c (x > Xc) to negative at t > t c (x < Xc)- As 
described in Ref. [44], the quantum evolution of the system around t c can indeed be solved exactly, as long as the 
non-linearities in the Higgs field (i.e. the self- interactions of tp), the interactions with the inflaton fluctuations and 
with the gauge fields are neglected. In such a case the Higgs behaves as a free scalar field with time dependent mass 
rritpit), whose quantum evolution can be solved in terms of Airy functions [44 . When t J> t c , the low momentum 
k < m Higgs modes grow exponentially fast due to the negative square mass, in a process known as "tachyonic 
preheating" [331 SI] [33J H2] . This in turn sources the production of the gauge field through its interactions with the 
Higgs. 



d X (t) 



d(mt) 



9\Xc 



(10) 
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The low momentum Higgs modes acquire large occupation numbers during tachyonic preheating and therefore 
become quasi-classical very fast. The transition from quantum to classical of the Higgs' long wavelength modes, takes 
place very rapidly indeed as compared with the time scale of the symmetry breaking process itself. This way, soon 
after the end of inflation but much earlier than the Higgs reaches the true vaccum, the modes within the tachyonic 
band have become fully classical, while all other modes remain in their quantum ground state. Due to the non- 
linear A0 4 interaction term, the non-tachyonic modes will eventually get populated thanks to the highly occupied 
low-momentum modes. Moreover, since the Higgs and the gauge helds are directly coupled, the gauge helds are also 
excited during and after tachyonic preheating. 

Our strategy here, following [¥TJ[H], will be therefore to introduce the system into a lattice in the precise moment 
in which the Higgs modes within the tachyonic band have became classical (slightly after the inflaton crossed the 
critical point) but the non-linearities and other interactions are yet negligible. In order to solve for the evolution of the 
system we will then use the lattice approach, consisting of replacing the quantum evolution by the classical one. This 
way, the quantum nature of the problem remains in the stochastic character of the initial conditions. This approach 
is, of course, only valid in the particular scenario under study because the quantum dynamics of the low momentum 
modes of tp drives the system into a regime of classical held behaviour, which ultimately justifies the use of the lattice. 
It is the subsequent non-linear classical behaviour of the symmetry breaking held that induces the growth of the 
inflaton and gauge helds, which then develop a non-trivial anisotropic stress-tensor which will source gravitational 
waves. This approximation has the clear advantage that it fully captures the non-linear and non-perturbative nature 
of the problem, and allow for the use of gauge fields in a relatively simple way. The full non-linear evolution of the 
system can then be studied using lattice techniques, discretizing both in space and time the classical equations of 
all fields @, © and @, see subsection |IV B| below. Note that other authors have also used the lattice formulation 
to study hybrid preheating in the presence of gauge helds, studying some of its phenomenological consequences like 
baryogenesis jJTJ 2H] or the generation of the primordial seed of cosmological magnetic fields [H| . Here, for 

the first time, we study the production of GW associated to the presence of such gauge fields. 

Finally, note that in general the model-dependent inflaton part of the potential Vmf (x) does not affect significantly 
the dynamics of preheating, except by setting the initial velocity with which the inflaton reaches the critical point. 
We will thus ignore VSnfl(x) and treat the initial velocity V c as a free parameter. 
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FIG. 1: Time evolution of the inflaton's mean (x)/Xc 
(blue) and of the Higgs' root mean squared /v (red) 
in a hybrid model with X — g 2 /2 — 0.125, e = 0.3 and 
Vc = 0.024. One can clearly see the growth of the Higgs 
towards the true vacuum, while the inflaton rolls down to 
the bottom of the potential. Once the Higgs is sufficiently 
close to the true vacuum, its self-interactions compensate 
the tachyonic mass and the field oscillates close to the 
VEV, while the inflaton oscillates around x = 0. 



FIG. 2: Time evolution of the covariant energy density of 
the Higgs (Diip(Diip)*) (black), the electric (EiEi) (red) 
and magnetic (BiBi) (blue) energy densities, and the in- 
flaton's gradient energy density (dixdiX) (green), all nor- 
malized to the total energy of the system. The parame- 
ters are the same as in Fig.[l] A fast growth of the energy 
components is experienced during tachyonic preheating, 
first in the Higgs field and then in the rest of the fields 
due to the their interactions with the Higgs. 



B. Initial Conditions 



The initial conditions of the fields follow the prescription from Refs. [UllUj. The Higgs modes ipk are solutions 
of the evolution equations obtained from ([2]), which can be rewritten for each component, as oscillators with time- 
dependent frequencies ipk + [k 2 — 2m 3 V c (t — t c )] Lp k = 0. If V c is not extremelly small and the couplings A, e and g are 
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not very big, see Ref. [H], there is always a time U greater than the critial time t c , but much shorter than the time 
scale in which the non-linearities are important, for which the Higgs modes within the tachyonic band have become 
classical, whereas those out of the tachyonic band can be set classically to zero. The amplitude of the tachyonic modes 
can then be found at the time t$ , as distributed according to a Gaussian random field of zero mean, which translates 
into a Rayleigh distribution 



P(\<p k \)d\<p k \dB k = m exp {-\<p k \ 2 /ai} d\ Vk \d9 k , 

7Ta k 



(11) 



for the modulus, with a uniform random phase 9 k € [0, 2ir}. The dispersion a\ is given by a\ = P(fc,^)/fc 3 , 
where P(k,U) = \f k \ 2 is the power spectrum of the initial Higgs quantum fluctuations in the background of the 
homogeneous inflaton [44 . In the classical limit, the conjugate momentum <^fc(r) is uniquely determined through 
4>k{U) = F(k,U)(pk(U), where F(k,r) = lm(f^(T)gk(r)), with f k and g k functions defined in Eqs. (43) and (44) 
of [H]. The rest of the fields (the gauge fields, the non-zero modes of the inflaton and the gravitational waves) are 
supposed to start from the vacuum, and therefore we will semi-classically set them to zero initially in the simulations. 
Their coupling to the Higgs modes will drive their evolution, giving rise to a rapid growth of the gauge and inflaton 
modes and GW subsequently. Their non-linear evolution will then be well described by the lattice simulations. 

One could think that, initially at time ij, since the vector fields are in vacuum (they have not been excited yet), 
one could put their amplitude and conjugate momenta to zero, once the classical regime in the Higgs sector has been 
established. However, the gauss constraint in Eq. ([6| imposes a relation between the gauge fields and the Higgs' 
components, which must be fulfilled at any time during the evolution. Therefore, given our initial conditions for the 
Higgs described above, we initialize the gauge field in such a way that Gauss constraint is satisfied initially, following 
[41] . In the temporal gauge Aq = 0, the Gauss constraint Q in Fourier space reads 

+ ik i A i (k,t i )=jo(k,t i ), (12) 

with jo = — 2e 2 Im Note that our initial conditions for the Higgs ensures that the initial charge vanishes, 

J d 3 xjo(ti) — 0. This is a necessary condition for Eq. (|9j) to be satisfied initially. We then set the initial conditions 
for the electric field in Fourier space, according to 

Ai(k,U) = -i^j (k,U) (13) 
Aiik = 0,ti) = (14) 



which is a particular solution of (12). Thus our initial conditions for the gauge field are such that initially the electric 



field is purely longitudinal, Aj(fj) oc fej, while the magnetic field vanishes, Aj(f,-) = 0. 

III. THEORETICAL PERSPECTIVES ON THE GW PRODUCTION IN SCALAR GAUGE THEORIES 

Our main purpose in this paper is to study the details of the production of a stochastic GWB during the reheating of 
the universe after hybrid inflation. Such GWB has indeed been extensively studied recently in global models of hybrid 
inflation, in the absence of gauge fields [28 31j. There the dynamics of the interacting fields generate a non-trivial 
anisotropic stress-tensor in the scalar fields, which source a significant background of GW. Here our aim is to embed 
the hybrid model described in Eq. into the gauge invariant framework of Eq. ([2]). The symmetry breaking field 
tp will then be coupled to the corresponding gauge fields associated to the gauge symmetry, and the dynamics of all 
the fields will be described by the set of differential coupled equations Q, ^ and ([9j. We will then study in detail 



numerically (see sections |IV| |V| and VI I the production of GW as coming form the new sources due to the presence 
of the gauge fields. Before going into the details of this model, however, we will first highlight some theoretical and 
model-independent results about the production of GW from gauge fields at preheating. 

In this section, therefore, we start laying down the basic formalism to study GW production in scalar gauge theories. 
Following Ref. |30j , we then apply this formalism to very simple, but relevant, wave- like sources which allow for an 
easy analytical treatment. We will review the argument why massless vector fields might enhance GW production 
during the turbulent evolution towards thermal equilibrium after preheating and we will extend this argument to 
gauge-invariant theories. We will then show that massless gauge fields are not produced at preheating in abelian 
scalar gauge theories. 

Gravitational waves on top of a (spatially flat) FLRW background correspond to (gauge invariant) linear pertur- 
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bations of the metric that are symmetric, transverse and traceless (TT), 8gij = a 2 (t) hij with dihij = ha = 5 . The 
perturbed Einstein equations describe the evolution of these tensor perturbations |52j as 

hij + ZHhn - W 2 hij = 16ttG Uj^ . (15) 

The source term is the TT part - verifying diUj? = nT T = - of the anisotropic stress, a 2 Ily = Tij —pgij where p is 
the background pressure. In the present case, both the scalar fields (Higgs and inflaton) and the vector fields (gauge 
bosons) contribute to the GW source. 

Defining — a and going to Fourier space, the GW equations ( 15 ) read 



+ [k 2 - J hy = 16tt G a 3 n^ T (16) 

where k 2 = k 2 is the square of the comoving wave-number and primes denote derivatives with respect to conformal 
time dr = dt/a. The TT part of a symmetric tensor in Fourier space is obtained by the projection 



Uj^{k) = O mm (k)U lm (k) (17) 
1 
2 



O ijtlm (k) = P a {k) P jm {k) - - Py(fc) P lm (k) , (18) 



where Pij(k) = 5ij — ki kj and k = k/k is the unit vector in the direction of k. 

For causal processes like preheating, most of the GW are produced well inside the Hubble radius, where the term 
in a" I a can be neglected in (16). The solution of (16 1 is then expressed in terms of a simple Green function 6 



h ij (k,r)= 16 ^ I ,/f sin [/,-(- -f I] ii''(t\ TT ' ( ' (k. r i (!«)) 



for initial conditions hij = h'^ = at r = Tj. 

The expansion of the universe will not play an important role for the discussion in this section and it is usually 
negligible during tachyonic preheating after hybrid inflation that we will study in the next sections. Therefore, from 
now on, we will work with a Minkowski background, a = 1. In that case, the TT part of the anisotropic stress for the 
abelian model (J2j) reads 

n^ T = [diXdjX + 2 Re [D iV > (D jV >)*] - B.B, - E t E 3 ] TT (20) 

where pure-trace terms have been removed by the TT projection. Non-abelian gauge fields would lead to similar 
contributions, with the corresponding modifications of the covariant derivative and the electric and magnetic fields. 



A. GW from wave- like sources 



In general, the calculation of GW production through ( 19 1 requires the knowledge of the time evolution of the 
source. It is interesting to study the very simple case where the source corresponds to a superposition of waves with 
wave-like dispersion relation and adiabatically evolving frequencies j3U]. This case covers to a good approximation 
different situations that occur in the context of preheating, in particular the stage of turbulent evolution towards 
thermal equilibrium 7 . Specifically, consider the following time evolution for the Fourier modes of some generic scalar 
field </>, 

0(fc,r) =a+(r)e lw *W r +a-(r)e- l ^« T with u 2 (k) = k 2 +ml (21) 



5 In the following, we will raise or lower indices of the metric perturbations with the spatial metric, which reduces at first order in h to 
the Kronecker 64 j, so hij = W = h}i and so on. 

6 This is also the exact solution for all the modes in a radiation-dominated universe where a" = 0. This is a good approximation in a 
wide class of preheating models where the equation of state quickly jump towards w = 1/3 at the beginning of preheating 53. • In a 
matter-dominated universe, the Green function for wavelengths of the order of, or larger than the Hubble radius is given by other Bessel 
functions. 

7 Of course, there are also several situations in the context of preheating where the time-evolution of the fields is not wave-like, such as 
the exponential amplification of fluctuations or the collision of non-linear bubble configurations. 



and similarly for the other fields. In the context of preheating, we deal generally with interacting waves where the 
interactions contribute to the effective mass m$. The cases in which the frequency is not constant but is nevertheless 
evolving adiabatically with time (ui -C w 2 ), can be treated in the same way, since there the WKB approximation gives 
<f>(k, t) oc exp (±i J T uj^^k, f ) df). 

In theories involving only scalar fields, the source terms for GW have the same form as the first term in the RHS 
of ( 20 ) . In Fourier space, the product of the spatial derivatives of the scalar field leads to the convolution 



n£ T (fc,-r) ocO ij!lm (k) 



d 3 p pi p m <f>(p, t) 4>(k - p, t) 



(22) 



where we have used Oij^ m k m = 0. Inserting (22) into (19), and using the time dependence (21), we get, whenever 



the coefficients a k (r) evolve adiabatically with time, 

hij(k,r) oc G e ±lkT O ijt i m {k) d 3 p Pip m oi^ af_ p I df exp i f I ±u}$(p) ± u^flfc - p\) T k 



(23) 



where we have decomposed the sine in ( 19 ) into imaginary exponentials. In the limit of large time r with respect 



to the frequencies, the time integrals above reduces to Dirac delta distributions, enforcing energy conservation for 
trilinear processes involving two "particles" of the field <fi and one graviton. For instance, we have 



u^ip) + uj<j,(\k -p\) = k . 



(24) 



for two ^-particles emitting a graviton (left panel of Fig. [3]). Not only the conservation of energy, but also the 
conservation of momentum is taken into account in ( 24 ) . Other signs in the phase inside the time integral in ( 23 ) 
correspond to other trilinear processes, such as the decay of a 0-particle into another ^-particle and a graviton. 
Energy and momentum conservations are possible only for massless ^-particles 8 and for k \\ p. However, for k \\ p, 
the TT projection brings the GW amplitude to PiPm = in rt23b. The reason for this is clearly helicity 



conservation, which forbids interactions between scalar waves and a graviton at leading order in the gravitational 
coupling constant G. Interactions involving several gravitons are possible beyond the linear order in hij, but they are 
highly suppressed by extra powers of the Newton constant G. 



<Kp) 



h H (k) 




* ( k-p) 

FIG. 3: Contribution of different source terms to GW production from wave-like fields. Left pannel: contribution associated to 
the source term in di(f>dj<j} for scalar fields, corresponding to the interaction between two scalar waves and a graviton (forbidden 



by helicity conservation). Middle pannel: contribution associated to the terms in ( 25 1, corresponding to the interaction between 



two vector waves and a graviton (allowed if the vector field is massless). Right pannel: contribution associated to the second 
term in (261, corresponding to the interaction between several scalar and vector waves and a graviton. 



On the other hand, the presence of massless vector fields (photons) may change this result, since interactions between 
two vector waves with helicity 1 and a graviton (middle pannel of Fig. [3]) are not forbidden by helicity conservation. 
Indeed, consider the two terms BiBj and EiEj in the source (20) of GW, and decompose the vector field Ai into a 
longitudinal part A\ = c^L and a transverse part Aj with diAj = 0. Proceding as above leads to energy-momentum 



That Eq. i |24[ ) implies = simply reflects the well-kown fact that a massless particle (here a graviton) cannot be emitted by two 
massive particles (here the </>-particles), since in the frame of the center of mass of the two incident particles the massless particle should 
be emitted with zero momentum. 
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conservation like in ( 24 ) , which again implies the momenta of the two vector particles to be parallel to the momentum 
of the graviton. As before, each term involving a spatial derivative di or dj then leads to a factor of fcj or kj in 
Fourier space, which vanishes when contracted with Oij t i m (k). However, the terms BiBj and EiEj in (20), include 
contributions as 



d h At dkAi and 



AJAJ 



(25) 



which do not involve di and dj. These terms can therefore lead to a non-zero GW amplitude from a wave-like source, 
through processes like the one illustrated in the middle pannel of Fig. [3] Note however that energy-momentum 
conservation for these processes still requires the vector field to be effectively massless. 

In addition, extra channels of GW production from wave-like fields arise in scalar gauge theories from the covariant 
derivative of the scalars. Indeed, consider the second term in the source (20) of GW. It involves contributions of the 
form 



edi<pAj<p* and e 2 A i A j \p\ 2 . (26) 

Since these source terms are cubic and quartic in the fields, GW production from these terms for wave-like fields 
correspond to interactions between three and four particles and a graviton, see e.g. the right pannel of Fig. [3] Clearly, 
these higher-order interactions are less restrictive as for instance they are now possible even for massive scalar and 
gauge fields. 

These qualitative arguments, although valid only for a very simple time-dependence of the source, illustrate that 
gauge fields can lead to new channels of GW production. In particular, they suggest that gauge fields can keep GW 
production active well after preheating, during the stage of turbulent evolution towards thermal equilibrium. Since 
this stage can last for a long period of time, this could significantly increase the amplitude of the resulting GW 
spectrum. As we will see, this does not happen for the hybrid model described by and (|2|). In the context of 
the hybrid abelian-Higgs models, the gauge field will have specific effects on GW production during preheating but 
the GW amplitude will eventually saturate during the later stages of the dynamics. Indeed, in such abelian model 
the gauge field acquires a mass through the Higgs mechanism during the dynamical spontaneous symmetry breaking. 
Thefore, as discussed above, processes like the one in the middel pannel of Fig. [3] are forbidden 9 . It is possible that 
massless gauge fields produced at preheating could still enhance GW production during the subsequent evolution 
towards thermal equilibrium. However, this does not happen in abelian scalar gauge theories, as we discuss next. 



B. No massless gauge fields from abelian scalar fields during preheating 

In the rest of this section, we argue that massless gauge fields are not produced during preheating in abelian scalar 
gauge theories, at least with canonical kinetic terms. More precisely, we show that in this case a would-be massless 
gauge field is either decoupled from the other scalar and gauge fields and therefore not produced 10 , or it receives an 
effective mass due to its interactions during preheating and the subsequent turbulent evolution 11 . Here and in the 
following, the term "mass" will be used in a very "loose" sense and should be understood as a contribution to the 
dispersion relation of the fields in the adiabatic regime (21), which is relevant for our discussion of GW production 
from wave-like sources in section III. A. We will consider a generic form of the scalar potential and an arbitrary number 
of U(l) symmetries and of scalar and gauge fields charged under these symmetries. However, the argument does not 
apply to non-abelian theories, like SU(2) x U(l), where massless photons may be produced at preheating |4"5] . 

Before considering the general case, consider again the system ^ of a single abelian gauge field A^ coupled to 
a single scalar field (p. A first possibility to obtain a massless gauge field would be that the VEV of (p vanishes, 
since otherwise the gauge field acquires a mass through the Higgs mechanism. However, in that case, A^ acquires an 
effective mass through its interaction with the fluctuations Sip of ip 

d u F liV +2e 2 \6(p\ 2 A li =2eIm[6( P *d l t5<p] . (27) 



9 On the other hand, processes like the one in the right pannel of Fig. [3] are in principle possible, although they are further suppressed by 
the gauge coupling constant. 

Here we do not consider fermions, which of course can also couple to gauge fields. Since the production of fermions is limited by Pauli 
blocking, they are not expected to lead to an abundant production of gauge fields at preheating. 

1 When the system equilibrates and cools down, such an effective mass will eventually become negligible, but this occurs at much later 
stages in the dynamics, much beyond our simulations' range. 
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In the second term of the LHS, the large and classical fluctuations of the scalar field behave as a classical VEV 
providing a mass to the gauge field by the Higgs mechanism. The fluctuations 5<p are produced whenever is 
produced, since both of them act as a source for the production of the other. Either the fluctuations of cp are 
responsible for the production of or the gauge field is produced by another mechanism and leads in turn to the 
production of Sip fluctuations through its interaction with (p. Thus it does not seem possible to produce a massless 
gauge field at preheating in this model. 

Another possibility is to consider several gauge fields so that a combination of them can remain massless like in the 
517(2) x £7(1) Standard Model. Consider therefore an abelian model with two gauge fields and A^\ coupled to 
a single scalar field ip with coupling constants e± and e 2 , D^tp = — ie\A^ — ieg-A^J (p. The system is invariant 
under the gauge transformation 

p _> e i«(*)+i/J(«) ^ , A W {x) ^A^{x) + -d ll a{x) , A^(x)^A^(x) + -d^(x)- (28) 

The equations of motion are 

d'F^J =2e k Im&*D^} (29) 

for k — 1,2 and where F^J = d^A^ — d v Ap . It is easy to see that the combination A^ = e 2 A^ — e\ A^ is 
massless. However, both the mass term and the source term vanish for A„, so it is a free field 

d v F^ = (30) 

where F^ v — d^A^ — d v A^ is the gauge field strength of the massless combination A^. Thus the massless gauge field 
is decoupled from the system and therefore not produced. This is similar to the photon of the SU(2) x U(l) Standard 
Model, which does not couple to the Higgs at tree level. However, it couples to the W and Z bosons because of the 
non-linear nature of non-abelian theories. In that case, Higgs fluctuations amplified at preheating may source the 
production of W and Z bosons, which in turn source the production of photons [32]. On the other hand, in the case 
of U(l) symmetries, interactions between gauge fields would require non-canonical kinetic terms in the lagrangian, 
like (F^F^) , or higher order terms like QF^ V W V or $ <pF^F^ p ', in order to avoid the tree-level decoupling of the 
photons. 

It is straightforward to extend the arguments above to an arbitrary number of scalar and gauge fields, see Ap- 
pendix |Xj We conclude thus that massless gauge fields are not significantly produced during preheating in abelian 
scalar gauge theories. As discussed in the previous subsection, we therefore do not expect a significant production 
of GW during the stage of turbulent interactions after preheating in such theories. The situation could be different, 
however, for non-abelian theories. 



IV. NUMERICAL CALCULATION OF GW WITH GAUGE FIELDS 



In this section, we present our numerical method to calculate GW production in scalar gauge theories on the lattice. 
The basic methods developed in [29H31] for scalar fields can be directly generalized to such theories. We briefly review 
these methods in the first subsection. On the other hand, in the presence of gauge fields, special care must be paid in 
the lattice calculation of GW in order to reproduce the continuum theory up to Oidx 1 ) and 0(dt 2 ) accuracy in the 



lattice spacing dx and timestcp dt. This is discussed in subsection IV B 



The numerical results presented in this paper were obtained with lattices of 128 3 and 256 3 points. We performed 
a number of checks to verify that the results are physical and not affected by lattice artefacts like insufficient IR or 
UV coverage or too large timestep. 



A. Numerical Method 



Since the expansion of the universe is negligible during preheating after hybrid inflation, we will work in a Minkowski 
background to simplify the notations. The extension to an expanding universe is straightforward |29H31j . In a 



Minkowski background, the equation ( 15 1 describing the evolution in time of the TT tensor perturbations representing 
GW reduces to 



d k d k hi 



167rGnT. T 



(31) 
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where the source term n^ T was given in ( 20 ) . In the lattice simulations, however, we can solve an alternative equation 
in position space, 



where the source term 12 



Hij - dkdkUij = 167rGII, 



Hij = dixdjX + 2 Re [D iV {D jV )*\ - B t B : - E t E 3 



(32) 



(33) 



is not Transverse- Traceless. We can thus solve eq. ([32]) together with the evolution equations of the scalar and gauge 
fields sourcing GW. The reason is the following: the TT part of the source is most easily calculated in Fourier space 
through the projection (17), which is non-local in position space. One could of course Fourier transform the anisotropic 



stress tensor at each timcstcp in order to calculate its TT part and then evolve the GW equation, but this would 
be highly time-consuming. However, since the equation (31), the TT projection (17) and the Fourier transform are 



all linear in h^, we can just solve for Eq. (32) with a non-TT source term, and apply the TT projection (17) on Uij 
in Fourier space only at those moments of time when we want to output the GW spectra. This is certainly a faster 
procedure since it does not require to take Fourier transforms at each timestep. This method was originally proposed 
in [29] . A modified version was developed in [31] based on the formalism of [30] . Another alternative method is to use 



the formal solution of ( 16 1 in terms of its Green function ( 19 ) to directly calculate the GW spectra from the source 

IE T I 



The results presented in the next section were o 



Dtained using the methods of both Refs. 



and 



and 



we checked that the method of [30] give the same results. We refer the reader to those references for more details. 
The GW propagate freely after their production and their energy density can be calculated by the spatial average 



Pgw 



(hij hij) 
32ttG 



(34) 



Performing an extra time average over a full period in order to eliminate the fast oscillations of the waves with time, 
the energy density is given by the sum of the kinetic and gradient terms 



1 



64nGV 



d 6 k 



(hij{k) hij(k)* + k 2 h^k) hij{ky 



(35) 



where V is the volume of the lattice box and h^ — Oijj m ui m is the TT component calculated in Fourier space. The 



spectrum of energy density in GW per logarithmic frequency interval at the time of production is then given by 

A- 3 



( dfgw] 
\<Mk), 



16GV 



(hi S {k) hijik)* + k 2 hij{k) hij(ky 



(36) 



where we have used the fact that the spectrum is statistically isotropic. 

Finally, the quantity of interest is the present-day spectrum of energy density in GW per logarithmic frequency 
interval divided by the critical density 



h 2 figw(/) 



1 dfgw\ 
Pc dlnf J Q 



(37) 



where the subscript "0" refers to today. Assuming a "standard" thermal history after reheating, the frequency and 
amplitude of the GW today are obtained from the spectrum at the time of production according to 



h 2 a, 



4 x 10 10 Hz 



9.3 x 10" 



1/4 
dp Pp 



i(l-3«i) 



I dp gw \ 
p dink J 



l-3lD 



(38) 
(39) 



where a p and p p = Aw 4 /4 are the scale factor and total energy density at the time of GW production, a* is the 
scale factor when the universe becomes radiation-dominated and id is the mean equation of state between these two 
moments of time, see e.g. |30j for details. Since tachyonic preheating is a very fast process, we assume that the 
universe becomes radiation-dominated in less than a Hubble time, so that the dependence on the scale factor in 
Eqs. p8p9] ) is negligible. 



Here we neglect again pure-trace terms in Uij because they will be removed by the TT projection. 
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B. Lattice Formulation 



Our lattice formulation for the evolution of the scalar and gauge fields is standard and described in full details in 
Appendix [B] The scalar fields are defined at the lattice points and the gauge field in the segments between lattice 
points. We start from a discretized version of the continuum action of the abelian-Higgs model ([!]) 



S 



-dtd 3 x ^2 
{x*} 



\ d+X d + ^x + {D+^Y + - + V( X , \<p\) 



(40) 



where the sum is over all the spacetime lattice points. We denote the lattice fields with a tilde to distinguish them from 
their continuum analog. The lattice expressions for the forward partial derivative d£, the forward gauge covariant 



derivative and the gauge field strength F^ are given in Eqs. (B2), (B3) and (B4) respectively. The action (40) 
is invariant under the lattice gauge transformation 



<p(x) -> e t&(x) <p(x) 



An 



An 



(41) 



which is a discretized version of the gauge transformation in the continuum. This ensures in particular that the 
discretized equations of motion derived form ( 40 ) lead to the conservation of the discretized version of Gauss con- 



straint |41) . These equations are derived in Appendix [Bj where we show that they reproduce the continuum equations 
up to 0{dx 2 ) and 0(dt 2 ) accuracy in the lattice spacing dx and timestep dt. 

On the other hand, this property would be lost by a naive discretization of the equations of motion (|3T|[32l for GW. 



The non-TT tensor perturbations Uij(x) are defined at the lattice points (as the scalar fields) and their discretized 
equation of motion is 



16ttGIL 



(42) 



Contrary to the equations of motion of the scalar and gauge fields, the lattice expressions d+ Z)+ and F^ v reproduce 



their continuum analog up to 0(dx, dt) only as dx, dt — ¥ 0. Therefore, the same would be true for the RHS of Eq. (42 ) 
if we use these lattice expressions to calculate the different source terms ( 33 ) of ITj . Instead of this, we have to 
construct new lattice expressions that reproduce their continuum analog up to 0(dx 2 , dt 2 ) accuracy and which lead 



to a stress-energy tensor that is still invariant under the gauge transformation ( 41 ) . The details of this procedure are 
derived in Appendix |B] and the final expression for 11^ is given in Eq. ( B23 1 . 



Note that the calculation of GW production with (D(dx 2 ) accuracy ensures a much better control on the UV part 
of the GW spectra. This is illustrated in Fig. [4j where we compare the GW spectra obtained by computing the GW 
source term with O(dx) and 0(dx 2 ) accuracy, for the same lattice and model parameters. The two spectra agree well 
in the IR, but the 0{dx) spectrum (in red) displays a larger (unphysical) growth in the UV. 
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FIG. 4: Example of GW spectra calculated with O(dx) 
(red) and 0(dx 2 ) (blue) accuracy (see Appendix |B|for de- 



tails), for the model (Tj2 1. The two spectra were obtained 
with the same model and lattice parameters and they 
were output at the same moment of time (mt — 100). 



FIG. 5: Evolution with time of the fraction of energy 
dens ity in GW during preheating, p gw /ptot, for the model 
( lpl with e/v^ = 6 (red), 0.5 (blue) and (black). The 
last case simply corresponds to the model without gauge 
field. The other parameters were v — 10 -3 Afpi, A = 
g 2 /2 = W' 4 and V c = 0.024. 
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Finally, we conclude this section by showing in Fig. [5] how the total energy density in GW is accumulated with 
time during preheating in the model ( lp l for typical values of the parameters. As anticipated in section III GW 



production eventually saturates at late times after symmetry breaking. Another observation that we can already make 
form Fig. [5] is that the total energy density in GW varies in a non-monotonic way with the ratio of the gauge coupling 
and the Higgs' self-coupling e/y/X. For e/y/X ~ 1, the gauge field leads to higher GW energy density than in the case 
with only scalar fields. This does not result from the mere addition of an extra field, but rather from an increase of 
the terms sourcing GW due to the dynamics of the coupled system of scalar and gauge fields. As we will see in section 
VI a crucial role in the dynamics of preheating in the model (Tj2| is played by cosmic string configurations of the 



fields in position space. When e/vA ~ 1, strings of the Higgs field and strings of the gauge field have the same width 



and lie on top of each other. In that case, the different source terms in (20 1 add to each others at the position of the 
strings. As we will now discuss, not only the amplitude of the GW signal is affected by the gauge field, but also its 
spectral properties. 

V. SCALAR, VECTOR AND GRAVITY WAVES SPECTRA 

In this section we present our results for the GW spectra produced by the system of scalar and gauge fields 
at preheating in the model ([T] [2]). In order to highlight the consequences of the presence of gauge fields, let us 
first quickly review the gross features of the GW spectra produced from preheating in models invovling only scalar 
fields. In that case, the spectra of the scalar fields amplified by preheating are usually strongly peaked around some 
typical momentum k*, which depends on the particular model considered and can usually be calculated analytically 
as a function of the parameters. The final GW spectrum depends essentially on this typical scale fc*, with a peak 
frequency and amplitude today that can be estimated according to 

/* « x 1010 Hz - ^ « 10-6 (f^y - ^ 

where H p and p p are the Hubble parameter and the total energy density at preheating when gravity waves are 
produced. The factor 10 -6 arises from the redshift of the GW radiation. In configuration space, i?* ~ 1/fc* corresponds 
to the typical size of the "bubble-shaped" fluctuations of the scalar fields amplified by preheating. Not surprisingly, 
similar estimates hold for the GW produced by bubble collisions in first order phase transitions, where in that case 
i? + is the typical size of the bubbles when they collide. 

As we will see, a major consequence of the presence of gauge fields is to introduce, in addition to /c*, new characteris- 
tic scales in the problem, which will be inherited by the final GW spectra. We will therefore first study the appearance 
of these scales in the spectra of the scalar and gauge fields themselves, before considering their consequences for the 
GW spectra. In the next section, we will see how these new scales arise from the dynamics of string-like spatial 
configurations of the scalar and gauge fields. Of course, at the practical level, the presence of different scales in the 
problem makes numerical simulations more difficult, since each scale has to be resolved efficiently in a single simula- 
tion. One can tune the parameters in such a way that the different scales coincide with each other, but important 
consequences of the gauge field can then be missed. It was therefore important for us to develop a lattice formulation 
accurate up to second order in the lattice spacing dx (see Appendix [B]) , as it allowed us to obtain reliable results for 
the UV behavior of the GW spectra while keeping a higher IR resolution as compared to a calculation accurate up to 
0(dx) only. Nevertheless, we will naturally be able to simulate only a restricted region of the parameter space. We 
will therefore study in some detail how our results vary with the model parameters, in order to extrapolate them to 
other regions of the parameter space. 

As discussed in section |Hj the model involves five independent parameters: the Higgs VEV v, its self-coupling A, 
its coupling to the inflaton g, the gauge coupling constant e and the initial velocity of the inflaton condensate at the 



critical point V c (10). The GW produced in this model without gauge field were first studied in [28-30 and more 
in detail, exploring the parameters space, in [51] . Without gauge field (e = 0), the GW spectra are already very 
sensitive to the remaining parameters. In general, however, very small values of the coupling constants A and/or g 2 
are required for these GW to have a sufficiently small frequency today to be observable. Neglecting the expansion of 
the universe (which is a good approximation for preheating after hybrid inflation, unless v is very high), the VEV v 
can be scaled out of the field equations and of the initial conditions by suitable redefinition of the fields and variables, 
so the dependence of physical quantities on v is known exactly. For the GW spectrum redshifted into present-day 
variables, one finds that the GW frequencies do not depend on v at all 13 , while the GW energy density scales as 



Despite the fact that different values of v lead to different values of the energy density p = A?j 4 /4 during inflation and preheating. 
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h 2 Vt gw oc Gv 2 . These scalings with v are preserved in the presence of gauge fields and we will simply take v = 10~ 3 Mpi 
throughout this section. Depending on the remaining parameters A, g and V c , three different dynamical regimes of 
GW production from preheating after hybrid inflation were identified in Ref. [3T|. In each regime, the scale and 
the resulting GW spectra vary in a very different way with the parameters A, g and V c . As we will see, the effects 
of the gauge field may also depend on which regime is considered. Our main interest here is on the consequences of 
the gauge field for GW production and, as could be easily expected, an important parameter in this respect is now 
the ratio e 2 /A. Note that, in a regime where a very small A is required for the GW to fall in an observable frequency 
range, the ratio e 2 /A may be huge. 

Let us first study the consequences of the gauge field for A ~ g 2 and a significant initial velocity V c . As far 
as the scalar sector is concerned, this is the easiest case to simulate as different dynamical scales are of the same 

1/3 

order of magnitude. In that case, the scalar fields are amplified with a typical momentum fc* ~ V c m, where 
m = vAf = is the mass of the Higgs' fluctuations (divided by y/2) around the minimum of the potential. In 

Fig. [gJ we show the spectra fc 3 \ip\ 2 and k 3 \B\ 2 of the Higgs and magnetic fields at late time, mt ~ 250 (when the field 
distributions have saturated and evolve very slowly with time), for A = g 2 /2 = 1(T 4 , V c = 0.024 and e/VX = 8. The 

1/3 

main observation is that the two spectra are peaked around well distinct scales, k ~ V c m for the Higgs and the 
vector mass fe~c» for the magnetic field. This is a rather unusual situation in the context of preheating with only 
scalar fields, where mode to mode interactions tend to smooth the differences in the spectra of the different fields. In 
the present case, the vector mass is the typical scale for the width and interactions of string-like spatial structures 
of the gauge field, see the next section. We also note a "UV bump" in the Higgs' spectrum, which is absent without 
gauge fields. As we will see, these features will be imprinted in the resulting GW spectrum. We found that the spectra 
of the gauge field are always peaked around its mass as long as it is more massive than the Higgs. This is illustrated 
in Fig. [7J where we plot the spectra of the magnetic energy density for different ratios of the vector and Higgs' masses, 
mA/m v = e/VzA. The spectra of the electric energy density behave in a similar way. This ratio plays a crucial role 
in the theory of cosmic strings (in the Abelian Higgs model), distinguishing between Type I (m^ > Tn v ) and Type II 
(jtia < intp) strings. For instance, the dynamics of multivortex solutions is governed by the fact that the interactions 
between vortices is attractive for tua > rriip and repulsive for niA < m v . Here we note that the spectrum of the gauge 
field tends to be peaked around the greatest of these scales, i.e. k ~ y\v if m v > ttia and k ~ ev if m v < ttia- Note 
however the "IR features" of the spectra for m v > itia- 




k/m k/m 



FIG. 6: Spectra fc 3 \<p\ 2 and k 3 \B\ 2 of the Higgs and mag- 
netic fields at mt = 250 for A = g 2 /2 = 10" 4 , V c = 0.024 
and e/vA = 8. Here the normalization of the amplitude 
of the two spectra is arbitrary and has been chosen only 
for convenience. The results were checked with N = 256 
simulations with km/m = 0.075, 0.1 and 0.15. 



FIG. 7: Spectrum of magnetic energy density per log- 
arithmic frequency interval divided by the total energy 
density, A oc fc 3 \B\ 2 , at mt = 250 for different 

values of e/vA. From left to right, e/vA = 0.2 (black), 
0.5 (blue), 2 (red), 4 (black, dashed) and 8 (green). The 
other parameters are the same as in Fig. [5] 



The presence of well-distinct characteristic scales for the scalar and gauge fields leads to specific signatures in the 
resulting GW spectra, as illustrated in Fig. [8] In order to compare the GW and matter fields' spectra, we will first 
consider the spectrum of energy density in GW per logarithmic frequency interval 

"!•(«- (sfe), < 44 » 

as a function of the wave-number k, both quantities being evaluated at the time of GW production (i.e. during 
preheating). Present-day redshifted spectra will be considered later on. We see in Fig. [8] the presence of three distinct 
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peaks in the GW spectrum for e 3> y/X: an IR peak around k ~ 0.25 to, a middle peak located around the Higgs 
mass, k ~ to, and a UV peak located around the vector mass k ~ et) (i.e. fc/m ~ e/y/X). We never encountered 
such features in models without gauge field, where the GW spectra are usually peaked around a single frequency (see 
e.g. the black spectrum in Fig. |9| even when different scales are present in the model. Contrary to the UV peak, the 
position and amplitude of the IR and middle peaks in Fig. [8] are independent of e/y/X, as long as e > y/X. We will 
see how they vary with the other parameters later on. The frequency of the IR peak tends to be smaller than in the 
case without gauge field, see Figs. § and § When e ~ y/X, the middle and UV peaks are superimposed. The resulting 
GW amplitude is greater than in the cases e 3> y/X and e <C VX, see Fig. [9J This is already the case for e/y/X = 6 in 
Fig. [i] (red spectrum), where the UV peak has a higher amplitude than for e/y/~\ = 8 (blue spectrum). Finally, for 
e <C y/X the GW spectrum becomes indistinguishable from the case without gauge field, see e.g. the black and green 
spectra in Fig. || for e/y/X = and 0.2 respectively. 
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FIG. 8: GW spectra |44| for e/y/X = 6 (red) and 8 
(blue) at mt = 250. The other parameters are the same 
as m Fig. [6] The results were checked with N = 256 
simulations with km/m = 0.075, 0.1 and 0.15. 



FIG. 9: Same as Fig.[i]for e/y/X = (black), 0.2 (green), 
0.5 (blue), 2 (red) and 4 (black, dashed). 



Up to now, we have considered the spectra of the scalar, vector and gravity waves at late times, when the dis- 
tributions have saturated and evolve very slowly with time. It is interesting to see how they build up with time, 
as for instance the different peaks in the GW spectra of Fig. [8] are actually produced at different moments of time. 
Furthermore, although the spectra of the magnetic or electric energy density are useful to single out the vector sector, 
the covariant gradient energy of the Higgs contributes to a greater fraction of the total energy density, see Fig. [2] 
and similarly the source of GW is dominated by the covariant gradient terms in (20). By analogy with |44|), we 



can consider the spectrum of the covariant gradient energy density of the Higgs, per logarithmic frequency interval, 
divided by the total energy density 



^ov(fc) = 



1 dp cov 
ptot dink 



(45) 



where p cov = (Ditp (Diip)*) is the covariant gradient energy density of the Higgs. 

The evolution with time of f?^ ov (fc) is shown in the left pannels of Fig. [To| together with the evolution of the 
spectrum of energy density in GW |44| in the right pannels. Each line in Fig. |10| corresponds to a different interval of 

1 /3 

time. From mt = 5 to 15.5 (first line), both spectra f2£ ov and f2| w are peaked around fc* ~ V c m and their amplitude 
increases exponentially with time as the tachyonic instability amplifies the Higgs' fluctuations. From mt = 15.5 to 
19 (second line), the peak of the spectra moves clearly towards UV, while their amplitude still increases. During this 
interval of time, the amplitude of the Higgs' fluctuations reach its vev in more and more regions of space and start to 
oscillate back in the potential. This leads to the collision of bubbles of the fluctuations and thiner and thiner regions 
of space where the Higgs is locally small and where it is energetically more favorable for the gauge field to be localized, 
as we will see in Section IvTl 

Next, from mt = 19 to 23 (third line in Fig. 10), the spectrum ^ ov moves back towards the IR, where two peaks 
form around fc* ~ 0.25 to and fc* ~ to. During this interval of time, the fluctuations of the Higgs have oscillated back 
in the potential and become small in more and more regions of space, where the gauge field tends to be localized. 
Indeed, we will see in the next section that the string-like configurations of the gauge field get significantly fatter 
during this interval of time, the increase of their width corresponding to the shift of the spectrum towards IR that we 
observe here. Meanwhile, the IR part of the GW spectrum starts to increase during this interval of time. Then from 
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mt = 23 to 26.5 (fourth line in Fig. 10 1, the two IR peaks in f2g ov disappear and the spectrum oscillates back towards 
UV. During this time interval, the amplitude of the Higgs' fluctuations goes back to its vev in more and more regions 
of space. We will see in the next section that the string-like configurations of the gauge field tend to fragment into 
smaller structures during this interval of time. We see here that the IR part of the GW spectrum grow significantly 
during this period. In fact, the IR peak of f2| w around fc* ~ 0.25 to reaches almost its final amplitude during this 
interval of time, while the UV part of the spectrum will still significantly increase. This IR peak is inherited from the 
peak of f2£ ov at the same location that formed during the previous interval of time and which has now disappeared. 

Finally, f2£ ov moves slightly back towards IR from mt — 26.5 to 29 (fifth line in Fig. 10 1 with a new peak forming 
around fc* ~ to, before going back to the UV from mt — 29 to 50 (sixth line in Fig. 10). During this interval of time, 
the GW spectrum increases significantly around fc* ~ to. From there on the spectrum of £l% ov has saturated and it 
is peaked around the vector mass, fc* ~ 6 m. The GW spectrum then slowly increases around this momentum to 
eventually reach its final form displayed in Fig. [8] (red curve) . 

To sum up, we see that the different peaks in the GW spectrum appear at different moments of time during the 
process of tachyonic preheating and symmetry breaking. They can be traced back to similar features in the spectrum 
of the covariant gradient energy density of the Higgs. However, whereas these features disappear in ^ ov , which 
becomes eventually peaked around the vector mass, the IR peaks remain in the GW spectrum, since GW decouple 
as soon as they are produced and their spectral shape remains unchanged since then. 

Let us now study how the GW spectra varies with the other parameters. Without gauge field, and still in the regime 
A < g 2 and significant initial velocity, the peak frequency of the GW spectrum at the time of production varies as 

k cx Vc^ 3 to and the peak amplitude as f2| w cx V c -2 ^ 3 v 2 , see Ref. [3T]. In Fig. [llj we show the GW spectra computed 

for V c = 0.024, g 2 = 2 A, e = 3\/A and two values of the Higgs' self-coupling: A = 10~ 2 and 10 -6 . The two spectra 
are almost on top of each others (in the units of k/m) so, as in the case without gauge field, the GW amplitude is 
independent of A while their frequency at the time of production varies as k cx to cx VA, as far as the ratios g 2 /X and 
e 2 /A are kept constant. The dynamics and the shape of the final GW spectrum are very sensitive to these ratios of 
the coupling constants but not to the absolute value of A 14 . 



In Fig. 
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we show the GW spectra computed for A = g /2 = 10 , e/v A = 3 and two values of the initial velocity 



of the inflaton at the critical point: V c = 0.024 and V c = 0.003. The frequency of the IR peak of the GW spectrum 
varies exactly as V^ 3 , as without gauge field. Its amplitude does not vary exactly as V c 2 ^ 3 , which even without 
gauge field is only approximate anyway, but nevertheless increases as V c decreases. Thus we see that the IR peak of 
the GW spectrum, which is present when e ^> V% varies with the parameters approximately in the same way as the 

GW spectrum produced without gauge field. For e ;> however, extra peaks are still present around the Higgs 
mass k ~ to and the vector mass fc~eu. These are on top of each other in Figs. [TT]and[T2j We see in Fig. [12] that 
the amplitude of at least one of these extra peaks seems to decrease when V c decreases, while the amplitude of the 
IR peak behaves in the opposite way. However, a quantitative estimate of the GW amplitude as a function of the 
intial velocity of the inflaton at the critical point would require to consider a larger range of values of V c . Since the 
frequency of the IR peak decreases with V c while the frequency of the two other peaks remains constant for e > yAj 
it becomes soon impossible to catch accurately the different scales in a single simulation. On the other hand, for 
e \/~X the gauge field becomes negligible and we recover the results of [31] . 

As discussed in |31) . the regime with A <J g 2 and a significant intial velocity V c that we have considered up to now 
is the easiest one to simulate but not the most interesting one from an observational perspective. Without gauge 
field, another regime of GW production in the model (Tib occurs for A <; g 2 and a small initial velocity V c , where the 
onset of preheating is driven by quantum diffusion of the fields around the critical point [31] , We will not study the 
consequences of the gauge field in this case here, but the main difference of this case compared to the previous one is 
in the onset of preheating, while the gauge field is amplified during the later stages. We thus expect the gauge field 
to have similar consequences as above. 

Finally, the third regime of GW production identified in [31] corresponds to the case g 2 <C A. In that case, the hybrid 
potential is much flatter in the inflaton direction than in the Higgs direction. As a result, the inflaton condensate 
oscillates several times around the minimum of the potential with relatively large amplitude, see Fig. |13| This leads 
to successive amplifications of the inflaton fluctuations by a combination of tachyonic instability and non-adiabatic 



resonance. In Fig. 14 we show the GW spectrum for X/g = 16 with (blue) and without (red) gauge field. As 



before, the gauge field leads to GW spectra with several peaks. We see here that the frequency of the IR peak can be 



Rescaling the coordinates as x — > rax and the fields as <j> — > <j>/v, the coupling constant A drops out of the equations of motion and 
enters only in the initial conditions for the amplitude of Higgs' fluctuations. These have very little consequences on the shape of the 
final GW spectrum, as illustrated in Fig. 1111 
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FIG. 10: Evolution with t ime of the spectrum of covariant gradient energy density of ip ([45} (left pannels) and of the spectrum 
of energy density in GW (441 (right pannels) for A = <? 2 /2 = 10~ 4 , V c = 0.024 and e/V A = 6. The spectra are shown between 
mt = 5 and 15.5 on the first line, between rat — 15.5 and 19 on the second line, between mt — 19 and 23 on the third line, 
between mt = 23 and 26.5 on the fourth line, between mt = 26.5 and 29 on the fifth line and between mt = 29 and 50 on the 
sixth line. In each case, the spectra are output every mt = 0.24 and red spectra correspond to earlier moments of time while 
blue spectra correspond to later moments of time (the colors go from red to yellow, green and blue as time evolves). 
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FIG. 11: GW spectra for A = 10" 2 (red) and A = 1(T 6 
(blue) at mt — 300. The other parameters are V c = 
0.024, g 2 = 2A and e = 3\/A- 



FIG. 12: GW spectra for V c = 0.024 (red) and V c = 
0.003 (blue) at mt — 250. The other parameters are A = 
g 2 /2 = 10" 4 and e/y/X = 3. The results were checked 
with N = 256 simulations with km/m — 0.05, 0.075 and 
0.1. 
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FIG. 13: Time evolution of the inflaton's mean (cj>)/ 'cj> c 
(blue) and of the Higgs' root mean squared {\<fi\)/v (red) 
for X/g 2 = 16, e/VX = 0.5, A = 10~ 4 and V c = 0.024. 
Notice the large oscilations of ((f>) as compared to Fig. [I] 



FIG. 14: GW spectra for X/g 2 = 16 and e/VX = (red, 
no gauge field) and e/yfX = 0.5 (blue) at mt = 400. 
The other parameters are A = 10~ 4 and V c = 0.024. 
The results were checked with N = 256 simulations with 
k IR /m = 0.03, 0.04 and 0.075. 



significantly lower than the frequency of the peak without gauge field. This prevents us from simulating a significant 
range of parameters, since both the IR and UV scales have to be included in the same simulations. Nevertheless, the 
results are consistent with a UV peak located around the vector mass, see Fig. 15 as was the case for A ~ g 2 . As 
before, we also expect the emergence of a middle peak located around the Higgs mass when e 3> VX, although we 
could not simulate this case for g 2 <C A. 



In any case, the main feature seen in Fig. 14 is the effect of the gauge field on the IR part of the GW spectrum, 
which can be of particular interest from the observational perspective. When X/g 2 increases, the frequency of the 
IR peak decreases while its amplitude increases, see Fig. [16} This is similar to the behavior of the GW spectrum 
without gauge field, which could be studied more accurately for a much wider range of values of X/g 2 in [3T]. It 
is thus possible that the gauge field merely enhances this behavior but that the frequency and amplitude of the IR 
peak remain relatively well described by the predictions of [31] without gauge field. On the other hand, it is perhaps 
not surprising to observe specific effects of the gauge field for g 2 <C A. As shown in Fig. 13 the inflaton condensate 
oscillates with relatively large amplitude around the minimum of the potential in that case. When it is far away from 
the minimum, the minimum in the Higgs direction is at \ip\ < v and the Higgs amplitude is small in relatively large 
regions of space. As we will see in the next section, these regions play a crucial role as magnetic strings are produced 
there. We therefore expect a different dynamics of the strings when g 2 -C A. A detailed study of this case is certainly 
interesting, but it lies beyond the scope of this paper. 

To conclude this section, let us summarize our results and relate them to the position space picture discussed in 
the next section. For e > VA, the final GW spectrum can be understood as made of three different peaks, which 
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FIG. 15: GW spectra for X/g 2 = 8 and e/^/X = 0.5 
(blue), 1 (black) and 2 (red) at rat = 350. The other 
parameters are V c = 0.024 and A = 10 -4 . The results 
were checked with N — 256 simulations with km/m = 
0.045 and 0.06. 
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FIG. 16: GW spectra for X/g 2 = 0.5 (black), 8 (red) 
and 16 (blue) at rat = 400. The other parameters are 
V c = 0.024, A = 10" 4 and e/VX = 0.5. The results were 
checked with N = 256 simulations with km/m = 0.03, 
0.04 and 0.075. 



are produced at different moments of time during the process of tachyonic preheating and symmetry breaking. A IR 
peak appears first, when bubbles of the Higgs start to collide and strings are formed in between the bubbles. The 
frequency of this IR peak tends to be smaller than the frequency of the peak of the GW spectrum produced without 
gauge field, but our results indicate that it varies in the same way with the model parameters. The frequency of the 
IR peak is set by the typical size of the bubbles of the Higgs field when they collide and by the correlation length 
of straight string segments at that time (these two scales are approximately equal). Next, a middle peak appears, 
whose frequency is set by the Higgs mass. This is the typical scale for the width and interactions of the Higgs' strings. 
Finally, a UV peak is formed when a significant fraction of energy has already been radiated away from the strings, 
see the next section. At that time, the spectrum of the gauge field is peaked around the vector mass, which is the 
typical scale for the width and interactions of the strings of the gauge field, and this scale sets the frequency of the 
UV peak of the GW spectrum. 



Once the GW are redshifted until today, their frequency is given by Eq. (38 1. We can then predict the present-day 
frequency of the three peaks as 

h £ f{9,\V c ) (IRpeak) 

/ 2 fa A 1 / 4 10 11 Hz (Middle peak) (46) 
/ 3 w A 1 / 4 10 11 Hz (UVpeak) 

where the frequency /i of the IR peak depends on the parameters g, A and V c and is of the order of, or smaller than, 
the frequency of the peak of the GW spectrum produced without gauge field. We can then use the predictions of [3T] 

A J vj 10 n Hz for g 2 £ A and V c £ 500 g 3 
f(g,\,V c )& { Ai. 9 10 n Hz for g 2 £ A and V c < 500 g 3 (47) 
A3 ^ 10 10 Hz for g 2 < A 

for the frequency of this peak. When e ~ vA, the middle and UV peaks merge into a single one with higher amplitude. 
For e <C vA, the gauge field becomes negligible and the results reduce to the case with only scalar fields, characterized 
by a GW spectrum with a single peak around / (g, A, V c ). Note that, depending on the parameters, the frequencies 
/i, f% and / 3 of the three peaks can differ by many orders of magnitude. It is of course not possible to simulate such 
cases on the lattice, but we expect the estimates ( |46"| to remain valid in such cases since they result directly from the 
presence of well-defined characteristic scales in the problem. 

Finally, we come to the amplitude of the GW spectrum. Without gauge field, or for e <C vA, we can use the 
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predictions of [21] 



10- 6 V C 2/3 (]^) 2 for g 2 > A and V c > 500 g 3 
i^ 8 (^) 2 for g 2 % A and V c £ 500 5 3 

(sf:) 2 g 2 « A 



(48) 



-5 A 



for the peak amplitude of the GW spectrum today. For e > vA, we could not vary the parameters over a sufficiently 
large range to study quantitatively how they affect the amplitude of each peak in the GW spectrum. However, our 
results indicate that the amplitude of the IR peak behaves roughly as in the case without gauge field and agrees 
within an order of magnitude with (48). The amplitude of the UV peaks should then be also well-described by these 



estimates if a significant fraction of the total energy density is indeed efficiently converted into small-scale structures 
of the Higgs and gauge fields. When e ~ vA, the middle and UV peaks are superimposed and their amplitude is 
slightly larger than in the case without gauge field. On the other hand, when e %/A, the amplitude of each peak is 
slightly smaller than in the case without gauge field. In any case, the amplitude of the three peaks is mostly sensitive 
to the VEV of the Higgs and the corresponding behavior is known exactly, f2 gw oc Gv 2 . 
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FIG. 17: Time evolution of the spatial distribution of the modulus of the Higgs field, |</?|(x, t) = \/ip\ + <p\, during the process 
of symmetry breaking. The images have been obtained with a N = 256 lattice simulation with an IR cut-off km = 0.15m, and 
parameters g 2 — 2A = 0.25, V c = 0.024 and e = 6 vA- From left to right, top to bottom, the snapshots correspond to mt = 
5.5, 11.0, 17.3 and 23.0. 



VI. SPATIAL CONFIGURATIONS 



The spectra of the anisotropic stresses of the matter fields and GW power spectra give us only a partial information 
on the evolution of the fields and the origin of the peaks in the spectrum. In order to understand the detailed dynamics 
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one has to use all the information available, and in particular, follow the spatial configurations in detail as a lunction 
of time, since then one can understand how specific features (like topological strings configurations) are formed and 
give rise to the observed peaks in the spectrum. Moreover, apart from both spatial images and power spectra, a very 
useful tool for this detailed understanding is the time evolution of histograms of both the Higgs and the magnetic 
fields' energy densities. These histograms allows us to identify the moments when the Higgs's oscillations make its 
vev reach zero and induce non-trivial windings at places where topological defects form. 

For this purpose we turn to the discussion of the production of GW in configuration space, describing the spatial 
distributions and correlations between the energy density of the scalar and vector sources, and that of the GW. We 
will see how the Higgs field forms bubbles that expand and collide and how the gauge field is excited during the 
symmetry breaking, forming elongated structures (tubes) in between the Higgs' bubbles. These string-like spatial 
configurations of the time-dependent Higgs and gauge fields follow from the dynamical equations of motion of the 
coupled system and will of course exhibit some differences with respect to the usual Nielsen-Olesen solution for static 
and infinite strings |54j . We will follow the formation and evolution of these strings during and after the symmetry 
breaking, but well before any scaling regime of string networks has been achieved and on length scales much smaller 
than the Hubble radius. 

For Nielsen-Olesen strings, the magnetic flux decays away from the core of the string with a typical length scale 
that is given by the inverse of the gauge field mass. However, during the process of symmetry breaking, the Higgs 
field provides an effective mass for the gauge field which oscillates in time as the Higgs oscillates around its VEV, 
with relatively large amplitude, see e.g. Fig. [T] Furthermore, the magnetic flux tends to be confined in the regions of 
space where the amplitude of the Higgs is small, i.e. between the bubbles of large Higgs amplitude, and the width of 
these regions is also strongly time-dependent as the Higgs bubbles oscillate and collide. Therefore, we will see that 
the width of the string-like configurations of the gauge field oscillate with time, becoming thicker when the amplitude 
of the Higgs is a maximum. The same occurs for the Higgs field whose effective mass is also time-dependent via its 
coupling to the oscillating inflaton field. At later times, we will see how magnetic energy density is shedded away 
from the core of the strings. Similar effects have been observed in field theory simulations of cosmological networks 
of cosmic strings in the Abelian-Higgs model, see e.g. [SS] and references therein. Such effects are usually suppressed 
by the ratio of the width of the string over the length of straight string segments, which is not small for the system 
under study. 




FIG. 18: Snapshots at time mt = 17 of the spatial distribution of the magnetic energy density B 2 (in units of m 4 ). The images 
have been obtained for a N = 256 lattice with an IR cut-off km = 0.1m, and parameters g 2 = 2A = 0.25, V c = 0.024 and 
e — 6 vA- In the left we see clearly how the string-like configuratinos of the magnetic fields are localized where the minima 
of the Higgs are, as described by the coloured transverse plane which plots the Higgs amplitude (normalized to the VEV) at 
that moment. On the right, the analagous figure where now the transverse coloured plane shows the phase of the Higgs, thus 
clearly demonstrating the existence of non-trivial winding in the Higgs around the positions of the the magnetic strings. 



A. Higgs bubbles 



In Fig 17 we show a sequence of snapshots of the spatial distribution of the modulus of the Higgs field for a model 
with couplings g 2 = 2A = 0.25, initial inflaton velocity V c = 0.024, and gauge coupling e = 6 vA. From left to right, 
top to bottom, the snapshots correspond to mt = 5.5, 11.0, 17.3 and 23.0, and capture the details of the symmetry 
breaking process towards the true VEV. We have choosen a simulation beginning at mt — 5, when the tachyonic 
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FIG. 19: Time evolution of the spatial distribution of the magnetic energy density B 2 (in units of m 4 ) along the process of 
the Higgs symmetry breaking. The images have been obtained with a N = 256 lattice simulation with an IR cut-off km = 0.1 
m, and parameters g 2 — 2A = 0.25, V c = 0.024 and e = 6 vA. From left to right, top to bottom, the snapshots correspond to 
mt = 5.5, 11.0, 17.3, 19.0, 21.0 and 23.0. At early times, before the Higgs bubbles percolate, the magnetic field is still very 
small and has not acquired yet the distinctive shape of topological string configurations. At times mt ~ 17— 19, the string-like 
spatial distributions of the magnetic energy density have finally developed, following the locus of points which corresponds to 
the intermediate regions between Higgs bubbles. The string-like distributions are most clearly seen at time mt — 19. Later, 
due to the time evolution of the gauge field's mass, the string segments fatten and start shedding away the magnetic field, see 
the main text. 
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modes of the Higgs are already well inside in the classical regime, as described by (11 1. Therefore, the distribution of 
the Higgs at mt = 5.5, i.e. sligthly after the initialization of the simulation, simply shows the bubble-like structures 
as developed in random positions, corresponding to the tail of high field values of the gaussian distribution ( 11 ) which 
describes the tachyonic modes. Since we are using a big lattice with N = 256 points per dimension, the number of 
bubbles we capture in a single box is statistically quite significant, and the resolution of their spatial profiles is also 
quite well captured, as can be clearly seen in the figures corresponding to mt = 5.5 and mt = 11. In the snapshot 
corresponding to mt = 11, the tachyonic instability of the excited modes has already led to the growth of the amplitude 
of the Higgs field towards the true vacuum, such that the Higgs field at the center of the bubbles at that time has 
reached already a 5% amplitude of the true VEV. The higher the random value of the Higgs was in a given location 
at the initial time, the faster the amplitude of the Higgs has grown in such location and the neightboring region. 
However, at time mt = 11, for the parameters choosen, the system has already entered into a regime in which the 
non-linearities due to the self-coupling for the Higgs are not negligible anymore. In the non-linear regime, the initially 
IR tachyonic modes are transferring power into the higher momenta modes out of the initial tachyonic band. As a 
result, in configuration space one can see that the amplitude of the Higss has grown everywhere in space, although 
the bubble structures are yet preserved. At mt = 17.25 we see that the Higgs has already reached the true VEV. 
Indeed, the Higgs at the centers of the bubbles has overpassed the VEV and reached a slightly greater value (as 
allowed by energy conservation). Due to the non-linearities, the bubbles are growing in that moment and are about 
to percolate. Since the size of the bubbles grows in time, the effective volume of the regions of lower amplitude in 
between the bubbles is naturally shrinking. At mt = 23, the bubbles have already percolated and one can clearly see 
in the intermediate regions between bubbles, that the Higgs amplitude is an order of magnitude lower than the true 
VEV. Those regions correspond to the locus of points in configuration space where a non-trivial winding has been 
developed, therefore leading to the formation of a spatial region where the Higgs amplitude is frustated to reach the 
true VEV. Of course, we are still far away from the stationary regime describing the evolution of topological defect 
networks. Rather, in the previous sequence of snapshots, we are looking at the dynamics of the symmetry breaking 
process itself, from the false to the true vaccuum, in a time scale in which the fields are still highly oscillatory and 
have not reached a scaling behaviour. 

The growth and collisions of the Higgs bubbles during symmetry breaking give rise to a significant anisotropic 
stress-tensor in the scalar fields, which drives the initial production of GW in this model. The specific correlation in 
configuration space between the locus of points where the gradients of the Higgs are maximum and the distribution 
of the GW energy density was already shown in [215] (section V of that paper) , see also [3T] . We will not reproduce 
again here such correlations and rather we will focus on the correlations between the GW energy density distribution 
and the energy density of the new source of GW considered in this paper, the gauge fields. 

Note that we will maintain the parameters choosen for Fig. [IT] fixed through the rest of this section, such that all 
the plots shown will represent some output from a simulation run with those values for the parameters. 



B. Magnetic string formation and evolution 

We already saw in momentum (Fourier) space how a new scale emerges in the spectrum of the gauge fields as the 
Higgs approches the true vaccuum, since then the gauge field aquires a mass through the Higgs mechanism. The 
covariant gradient energy of the Higgs is minimized for |41j 

e 

where f2 = 0/|0| is an element of U(l). This induces the magnetic field to concentrate its energy density within those 
regions in which the Higgs amplitude is smaller and phase gradients are larger. Thus, during the symmetry breaking 
of the Higgs, magnetic flux tubes will be concentrated in the regions between the Higgs bubbles. This can clearly 



be seen in the left pannel of Fig. 18 where we plot a high value isosurface of the magnetic energy density, showing 
this way how the three dimensional configurations of the magnetic field forms flux tubes (strings) . Together with the 
strings, we also show a two dimensional plane orthogonal to those long flux tubes, showing how the strings appear 
in those places where the Higgs is minimum. This correlation between Higgs' zeros and magnetic field strings is a 
universal feature of the Abelian-Higgs model, and is related to the Nielsen-Olesen vortices (and strings) predicted in 
the model many decades ago [54] . 

Note that in the Abelian-Higgs model, the magnetic flux along the strings is quantized and is related to the Higgs 
winding number, 

<I> , ; = / d 2 x B ■ z = (j) dx ■ A = — — a , 
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where n is an integer called the winding number. This topological number is conserved along the evolution unless 
there are over-the- barrier transitions during preheating, see Ref. |41j . In the right pannel of Fig. |18| we have shown 
the magnetic strings (the magnetic energy density) from a perspective in which one can see the (color coded) variation 
of the phase of the Higgs from to 27r, as painted in a transverse plane to the string in the center of the box. Such 
plane correspond to a bidimcnsional cut of the three dimensional distribution of the values of the Higgs' phase between 
and 2ir. Around the place where the central string segment touches the plane, one can clearly see that the Higgs 
phase winds non-trivially around the string. The plane is also crossed by another string segment close to one of the 
walls of the box, and again one can see the non-trivial structure of the winding around such string. Besides, the 
two-dimensional distribution of the Higgs phase in the choosen plane, also shows very clearly the locus of points (lines 
within the plane) where the Higgs phase jumps from 2-7T to 0. If we interpret the plane as a Riemann surface, the 
curve lines where the Higgs phase jumps would then be the edges of the Riemann surfaces. 

In Figure [19] we show the time evolution of the spatial distribution of the magnetic energy density, from time 
mt = 5.5 till mt — 23. As explained above, the magnetic fluxes will tend to concentrate in those regions in which 



the Higgs amplitude (gradients) will be minimum (maximum). Thus, in the first two plots of Fig. 19 we see that 
during the initial times of symmetry breaking (when the Higgs bubbles have not yet percolated), the distribution of 
the magnetic field forms inhomogeneous lumps where the magnetic energy density is maximum in the spatial regions 
in between the initial nucleated Higss bubbles. However, the amplitude of the magnetic field is still very small to 
compete with the gradients of the Higgs, so the GW production is driven initially only by the Higgs inhomogeneities. 
When the Higgs bubbles percolate, the magnetic field is finally excited significatively and its amplitude grows by 
several orders of magnitude. At the same time, the regions between the percolating bubbles shrink, forming elongated 
tubes, i.e. topological defects, which correspondingly induce new spatial configurations in the distribution of the 
magnetic field. In particular, at times mt = 17 and mt = 19, see the third and fourth plots of Fig. [19] the magnetic 
field is compressed into thin tubes, located precisely in the locus of points where the Higss was prevented to reach 
the true VEV due to the development of a topological winding around the tube. Since the effective mass of the gauge 
field oscillates dynamically according to the amplitude of the Higgs around the VEV, the string-configurations of the 
gauge field do not have a constant width, but rather a time-dependent one, see fifth plot of Fig. [19] Indeed, the 
magnetic field string-configurations get thicker and eventually break into concentric layers which are shedded away, 



see the sixth plot of Fig 19 Nevertheless, the Higgs winding remains, following the Higgs' zeros at the cores of the 
strings. At time mt = 23 we can even see in the figure a thin string at the core of a thicker one. Later, we will see 
that these features are also inherited by the spatial distribution of the energy density of GW. Since the magnetic field 
energy is pushed away from the core of the strings, one expects that, at late times, there will be three components 
of the magnetic field, one associated with the core of the strings, another one which has "evaporated" from the core 
and now occupies the whole box in the form of very small-scale structures, plus a diffuse radiation component [42] . 
This is indeed what we see in the first two plots of Fig[22] Again these features will also be inherited by the GW 
distribution of energy, see the last two plots of Fig. [2T] 

C. Histograms of the Higgs and magnetic field energy 

Histograms give yet another perspective on the dynamics of the Abelian- Higgs model at preheating. Power spectra 
gave us an idea of the energy distribution as a function of scale, and allowed us to pick up certain typical scales of 
the problem. Then we could look at the spatial distributions/configurations to search for specific features and find 
those given scales, like the string width and length. However, neither the power spectra nor the spatial distributions, 
tell us how common those features are. For that we have to look at the histograms of the values of the fields, and see 
how these distributions change with time. For instance, since magnetic strings seem to the localized around the zeros 
of the Higgs, if we can follow the time evolution of the Higgs histograms, we can see how often the Higgs field has a 
significant fraction of lattice points with its vev close to zero. At those points one expects new windings to appear, 



and indeed this is what it seems to occur, see left pannel in Fig. 20 

One can also correlate the time evolution (oscillations) of the Higgs histograms with the production of magnetic 
field energy, associated with both the strings and the diffuse component coming from their evaporation. We show the 



time evolution of these histograms in the right pannel of Fig. 20 



D. Spatial distributions of GW 

Finally, let us consider the spatial configurations of GW. The expression for the GW energy density was given in 
(34 1 and here we will look at the spatial distribution of hijhij. We will show in a sequence of snapthots the time 



evolution of |/i«(x, t)| 2 in space, similarly to what we did with the Higgs and magnetic fields previously. Note that 
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FIG. 20: Evolution in time of both the histograms of the Higgs field normalized to its VEV (left) and the magnetic energy 
B 2 normalized to m 4 (right). The first row corresponds to the initial times, from mt — 5.05 to mt — 19. The second row 
coresponds to times from mt — 19 to mt = 32. The third row corresponds to times from mt — 32 to mt — 50. It can be 
clearly distinguished that the Higgs moves very fast towards the true VEV in the intial stages of hybrid preheating and later 
oscillates close to VEV (see the left tails of the Higgs histograms). At any moment, even when the Higgs is oscillating in the 
broken phase with small amplitude compared to the VEV, there remains a significant fraction of points in the lattice where 
the amplitude of the Higgs is much smaller than the VEV, corresponding to the spatial regions located in between the Higgs 
bubbles. 
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FIG. 21: Time evolution of the spatial distribution of GW \hij (x, t)| 2 (in arbitrary units) along the process of the symmetry 
breaking of the Higgs. Note that we use the same arbitrary units for all the plots, so the relative amplitude between one 
snapshot and another tells about the physical growth of the GW energy density. The images have been obtained with a 
N = 256 lattice simulation with IR cut-off km = 0.1m and parameters g 2 — 2\ — 0.25, V c = 0.024 and e = 6 \/\. From left to 
right, top to bottom, the snapshots correspond to mt = 5.5, 11.0, 17.3, 19.0, 21.0 and 23.0. At early times, before the Higgs 
bubbles percolate, the GW energy density is still very small and distributed in lumps over the lattice, with maximum values 
in the regions where the gradients of the Higgs are maximum. At times mt ~ 17 — 19, however, the string-like configurations 
of the gauge and Higgs fields induce similar string-like distributions of GW. The tubes of highest energy in GW can be seen 
most clearly at times mt = 17 — 19. Later, due to the time evolution of the strings, i.e. their fattening and shedding away of 
small-scale structures, the distribution of GW seems also to follow a similar pattern. Particularly noticeable here is the figure 
corresponding to time mt = 21, where inside one of the concentric tubes one can see another string-like configuraion in the 
core of the tube, very similar to what we observed in the magnetic energy density. Finally note that the orientation of the box 
have been choosen such that one can clearly see as best as possible some of the features developed by the sp atia l distribution 
of GW. Consequently, there is no specific correlation between the particular magnetic strings shown in figure ( 19 1 and the ones 
shown here, since the boxes are being observed from very different orientations. 
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FIG. 22: Here we show the magnetic field energy density B 2 in units of m 4 (top) and the GW spatial distribution h 2 j in 
arbitrary units (bottom), for times mt = 29 (left) and mt — 35 (right). Note that there are in fact three components in the 
magnetic field spatial distribution, the very core of the strings (barely seen here), as well as lumps of a small-scale structure 
which has been shedded away from the initial string-like configuration, plus a diffuse background (in green in the figures) which 
we interpret as radiation. The spatial distribution of GW follows a similar pattern. 



we will plot |/ijj(x, t)| in arbitrary units, since we only want to hint the spatial features developed in the distribution 
of GW over the lattice, to show that these features follow precisely the string-like configurations of the source. As 
expected, the spatial distribution of GW will follow that of the dominant source term at each time. 

In Fig. 21 a series of six snapshots of \hij (x, t)\ 2 are shown sequentially from mt = 5.5 to mt — 23. In the first two 
figures, corresponding to times mt — 5.5 and mt = 11, the distribution of GW follows the spatial configurations of the 
Higgs bubbles during the initial tachyonic stage. The correlation between the Higgs features and the GW distribution 
was studied in detail in [55], here we simply want to emphasize that initially the spatial configurations of GW are 
distributed as lumps over the lattice. However, in the middle figures of Fig. |21| corresponding to times mt = 17 and 
mt = 19, the GW distribution has begun to concentrate within elongated regions which coincide precisely with those 
positions in space in which the Higgs and gauge fields have formed string-like configurations. The spatial distribution 
of GW is clearly concentrated around these strings. In the last two snapshots of the sequence of Fig. [21] corresponding 
to times mt = 21 and mt = 23, we see that the string-like configurations of GW also fatten and break into small 
structures shedded away in the lattice, in full analogy with the behaviour of the magnetic energy density that we 
observed before. At later times, we show in Fig. [22] how both the magnetic field and the GW are distributed all over 
the lattice, in the form of very small-scale structures. This behavior of the spatial distribution of GW, closely tracking 
the behavior of the abelian-Higgs strings as they are formed, evolve and later fragment into small-scale structures, is 
in perfect agreement with the successive appearence of the different peaks in the GW spectrum that we observed in 
Section V. 
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VII. DISCUSSION AND PERSPECTIVES 

Gravitational waves are a robust prediction of general relativity. There is indirect evidence of their existence 
from inspiraling binary pulsars, although no single direct detection has been claimed. A stochastic background of 
GW may soon be discovered, either directly with laser interferometer antennas or indirectly through the pattern of 
polarization anisotropies they induce in the cosmic microwave background. Such a detection would open a completely 
new and unexplored window into the Early Universe, possibly as rich as that which has been recently revealed in 
the CMB. There are many sources of GW that can generate a stochastic backgrounds and thus it is necessary to 
characterize those backgrounds with as much detail as possible. Apart from known astrophysical point-like sources 
beyond the confusion limit (where we cannot resolve them), there are also predictions for GWB from cosmic defects 
and hypothetical strongly first order phase transitions in the Early Universe. Moreover, Cosmological Inflation makes 
a robust prediction of a stochastic GWB produced during the quasi-exponential expansion of the Universe, with very 
specific spectral signatures: a Gaussian, almost scale-invariant spectrum with an amplitude directly related to the 
energy scale of inflation. If the scale is high enough (close to the Grand Unification Theories' scale) then these GW 
will also leave an imprint in the (curl) polarization anisotropies of the CMB. Unfortunately, this GWB is still too 
weak to be discovered with the near-future GW interferometric antennas, although it does cover a sufficiently broad 
frequency range to be detectable by future GW Observatories (GWO) like BBO or DECIGO. 

Furthermore, a robust prediction of inflation is that it must have ended, converting the vacuum energy responsible 
for the tremendous expansion into the matter and radiation we observe today. Such a process, known as reheating, 
is typically very violent and very inhomogeneous, with large density waves moving at relativistic speeds and colliding 
among each other, thus converting a large fraction of their gradient energy into gravitational waves. In some cases, 
the conversion is so sudden and violent that a significant fraction of the total energy that goes into radiation ends 
in a stochastic background of GW, which could be detected in the future. The energy spectrum of such a GWB is 
very non-thermal and far from scale invariant, but actually peaked at a frequency which is related to the typical mass 
scale responsible for the end of inflation (either the mass of the inflaton or that of the field that triggers the end of 
inflation, like in hybrid models), which could be orders of magnitude smaller than the Hubble scale at the end of 
inflation. However, if the energy scale of inflation is large (of order the GUT scale) then this stochastic GWB will 
be peaked at GHz frequencies, far from the present sensitivity of GW interferometers. Nowadays, our only chance 
of detecting the GWB from reheating is to consider the low-scale models of inflation - like hybrid models - with the 
appropriate parameters to convert a large fraction of the initial vacuum energy into GW. The analyses done so far 
have considered only scalar fields whose gradient energies source the anisotropic stresses needed for GW production. 
However, vector fields (gauge or not gauge) are expected to be an even better source of GW, due to their anisotropic 
curl components, so that preheating scenarios with gauge fields may have a larger contribution to the GWB than 
scalar models. In fact, previous studies of gauge fields at preheating, in the context of Electroweak Baryogenesis and 
in the generation of Primordial Magnetic Fields, have identified the formation of long-wave semiclassical gauge field 
configurations like sphalerons and helical magnetic flux tubes which evolve with time very anisotropically and could 
contribute significantly to the production of GW. 

In this paper, we have developed a formalism to calculate the production of GW by coupled systems of scalar and 
gauge fields on the lattice. The numerical method that we have constructed can be applied to different sources of 
GW where out-of-equilibrium gauge fields play an important role, such as thermal phase transitions, cosmological 
networks of local defects or non-perturbative decays of scalar condensates into gauge fields. We have studied in detail 
the dynamics and the production of GW during preheating after hybrid inflation, in the context of abelian-Higgs 
models that go through dynamical symmetry breaking triggered by the expectation value of the inflaton field. As 
the inflaton is driven slowly (as opposed to a quench) below the critical value, the mass squared of the abelian Higgs 
field becomes negative and drives the spinodal growth of long-wave modes of the Higgs. Since the Higgs is charged, 
its rapid growth induces a corresponding growth of gauge field configurations. At the end of inflation there are 
no temperature fluctuations that can induce over-the-barrier transitions. However, long-wave quantum fluctuations 
become semi-classical and act as a stochastic force that allow transitions over the false vacuum and thus induce 
(locally) the generation of a topological winding number of the Higgs field. After symmetry breaking, there is not 
enough energy to unwind the Higgs phase, leaving behind a Nielsen- Olesen string. Such cosmic string configurations 
can be seen explicitly in our spatial distributions of both Higgs and gauge fields. They play a crucial role in the 
production of GW at preheating and we observe that the spatial distribution of GW is indeed concentrated around 
the strings. Those strings will eventually decay (we see that they become wider and disperse away their energy density 
in the form of small-scale structures of the fields, although the winding phase around the core of the string remains, 
since it is topologically stable), which eventually shuts-off the GW production. 

The complicated dynamics occurring at preheating in this abelian Higgs-inflaton model has been studied using 
both power spectra analyses, as well as the fields' distributions in configuration space, together with histograms of 
the fields' values, as a function of time, in order to correlate the different features observed and their evolution. The 



29 



picture that arises is the following. At the end of inflation the tachyonic growth of the Higgs Gaussian random field 
creates an inhomogeneous distribution of fields characterized by "bubbles" of Higgs energy density that expand and 
collide. The gauge field concentrates at the valleys between the bubbles, where the Higgs has low values, forming long 
flux tubes of magnetic energy density. The dynamics of the bubles when they expand and collide leads to regions 
in space where the Higgs field reaches the false vacuum and there are over-the-barrier transitions, with topological 
windings associated with them. These Nielsen-Olesen vortices are connected with each other in a cosmic string which 
runs along the core of the magnetic flux tubes. There are strings that encompass the whole simulation box and even 
beyond, thanks to periodic boundary conditions. We observe this process both in configuration space and with the 
histograms of Higgs vevs. 

We have followed the dynamics of the strings during and after the symmetry breaking, although still on time scales 
shorter than the Hubble time and on length scales smaller than the Hubble radius. Once the strings are formed, they 
evolve by increasing their size and shedding away layers of magnetic energy density. At the cores of the strings there 
always remain a thin magnetic flux line but the energy seems to pour away from the strings in the form of waves 
concentric with the string. Nevertheless, we observe (in a transverse plane to the string) that at the core of the string 
there remains a conserved winding number of the Higgs. We have followed this winding number up to long times 
and we confirm that it is still there, in spite of the fact that the magnetic flux tube is so dilute that we cannot see 
it coherently: it seems to have "evaporated" . What remains is a diffuse background of small-scale structures of the 
Higgs and gauge fields permeating the whole box, together with the remnants of the strings. 

The formation, evolution and fragmentation of the strings are accompanied by a significant production of grav- 
itational waves which inherit specific features from the string dynamics. In position space, we observe how the 
distribution of GW follows very closely the evolution of the strings, being first concentrated around the straight seg- 
ments of strings, then fattening as the strings become wider and finally being dispersed over the lattice as the strings 
emit small-scale structures of the fields. In Fourier space, this dynamics is encoded into the successive appearance 
of very distinct peaks in the GW spectra. The position of each peak is directly related to the physical scales in the 
problem: the Higgs mass, which governs the width and interactions of Higgs field's strings, the gauge field mass, 
which governs the width and interactions of gauge field's strings, and the typical momentum amplified by tachyonic 
preheating, which determines the characteristic size of the bubbles when they collide and the correlation length of the 
straight segments of strings. The former two determine the peaks in the high momentum (UV) range of the spectrum, 
while the latter corresponds to the long- wave (IR) peak. The IR peak appears first, when the bubbles collide and 
the strings are formed, while the UV peaks are formed later on, when the strings evolve and decay into small-scale 
structures of the Higgs and gauge fields. When the different scales are close to each other, the different peaks are 
superimposed and the amplitude in GW increases. When the gauge coupling constant is significantly smaller than 
the Higgs' self-coupling, the results reduce to the GW spectra produced without gauge field, characterized by a single 
peak. 

We have calculated the GW spectra produced in this abelian Higgs model of preheating after hybrid inflation with 
state-of-the-art simulations, although still limited in spatial resolution and box sizes. In order to probe reliably the 
different scales in the problem in each simulation, we developed a lattice calculation of GW production with gauge 
fields that is accurate up to second order in the lattice spacing. Our numerical results for the GW spectra today 

). The present-day frequency and amplitude of these GW are very sensitive to the 
model parameters and the frequency of the different peaks may differ by many orders of magnitude, as illustrated 
in Fig. [23} As in the same model with only scalar fields, very small coupling constants are still neccessary for these 
GW to fall into a frequency range that is accessible by interferometric experiments. Whether this is natural or not 
depends on the underlying theory for inflation and particle-physics models of hybrid inflation with such small coupling 
constants have indeed been already proposed in the literature, see (31] and references therein. We also observed that 
the frequency of the IR peak in the GW spectrum can be smaller than the peak frequency produced in the same 
model with only scalar fields, so the gauge field may enlarge the regions of the parameter space that may lead to an 
observable signal. More generally, there are many other models of inflation and preheating where gauge fields may 
play an important role and wich may lead to GW that could be observed in the future. 

After preheating the system enters into a turbulent regime where, at least in the abelian Higgs model that we 
considered, gravitational waves are no longer produced and the GW energy density saturates. We expect this result 
to be rather generic in abelian scalar gauge theories, because in that case the gauge fields that are produced at 
preheating acquire a mass, either directly through the Higgs mechanism or due to their interactions with scalar fields' 
fluctuations. It would be interesting to study models with other gauge groups, like SU(2) x U(l), where gauge fields 
remain effectively massless after symmetry breaking. This may happen for instance during preheating after hybrid 
inflation close to the electroweak scale (possibly a secondary stage of inflation, not necessarily related to the CMB 
anisotropies, only responsible for reheating the universe), where the photon spectra may exhibit inverse cascade during 
the turbulent evolution towards thermal equilibrium |42j . This could significantly lower the typical frequency of the 
resulting GW today and relax the conditions on the parameters for these GW to be observable. The details of the GW 



are well described by Eqs. (46|48 
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FIG. 23: The predicted stochastic background of GW from preheating in the abelian-Higgs model for two different sets of 
parameters, A ~ 1CT 3 , g ~ 1CT 8 , e ~ 0.1, v ~ 10 11 GeV (red curve) and g ~ y/\ ~ 1CT 6 , e ~ 10 -4 , t> ~ 10 13 GeV, and negligible 
initial velocity (blue curve), together with the expected sensitivity from future GWO like Advanced LIGO/VIRGO, LISA, ET, 
BBO and DECIGO. N ote that in order to reach GWO sensitivity we had to extrapolate the position of the IR peaks using 
the expressions in Eq. ( |47[ ), and make an educated guess for the shape of the spectra between the peaks (dashed lines). Also 
plotted are the expected GWB from a global phase transition and that from an inflationary model with tensor to scalar ratio 
r = 0.1. 

spectra produced from preheating should also be rather sensitive to the particular gauge group under consideration 
because this determines the nature of the defects that can be formed. The defects do not have to be stable since these 
GW are produced when they are being formed. 

An intriguing possibility is the following. Given that preheating is so extremely inhomogeneous, and since these 
inhomogeneities get imprinted in the gravitational wave background, which immediately decouples from the plasma, 
one may envision a stage of technological development in the not so far future in which GWO with sufficient angular 
resolution may resolve the structures that gave rise to the GWB right at the moment when the Universe reheated. In 
the usual preheating scenario at high energy scales with only scalar fields, the physical structures will have today a size 
that is completely undetectable when projected over the sky, and thus the GWB will look essentially homogeneous 
from Earth. However, if gauge fields with long string-like configurations (of horizon size and possibly even with 
super horizon correlations) were behind the generation of the GWB, then one could expect to see inhomogeneities 
in the angular distribution of those gravitational waves. In particular, an array of GWO could detect the string- 
like anisotropies in the GWB across the sky. At the moment, the angular resolution of LIGO is not better than a 
degree projected in the sky. However, in the future one could resolve much finer structures in the GWB thanks to a 
dense network of ground-based laser interferometers, as proposed e.g. in Ref. [56] . Thinking ahead of our times, it 
may not be unrealistic to imagine that in the not so far future the GWB will be mastered with sufficient detail to 
resolve the anisotropies in this elusive background and thus recover vital information about the physics responsible 
for the violent conversion of energy from inflation to a radiation and matter dominated epoch (the Big Bang of the 
Old Theory). No other probe can give us so much information, since GW decouple immediately upon production 
and thus retain the spatial and energy distributions of the sources that produced them. We can compare with the 
Cosmic Microwave Background, which gives us detailed information about the epoch of photon decoupling thanks to 
the exquisite measurements of the angular correlation of both temperature and polarization anisotropies. The CMB 
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provides a snapshot of how the Universe was like 380,000 years after the Big Bang. On the other hand, the GWB 
would open a window into the physics of the Big Bang itself, allowing us to infer from its detailed features whether it 
was as violent and inhomogeneous as we predict, to determine what kind of fields were present and whether the rich 
phenomenology that we associate with preheating (topological defects, baryogenesis and/or leptogenesis, primordial 
magnetic seed creation, non-thermal production of dark matter, etc.) was actually realized in nature. 
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Appendix A: No Massless Gauge Fields in Abelian Scalar Gauge Theories During Preheating 



In this Appendix, we extend to an arbitrary number of fields the argument of Section III that massless gauge fields 
are not produced during preheating in abelian scalar gauge theories. 

Consider an abelian theory U(l) K , with K gauge fields coupled to TV complex scalars ip n 

K N 



- C = \ E F $ F(k) + E D *Vn (Al) 
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D^n = dpifn - i ^ e kn AW tp n (A2) 
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and Fj£) — d„Au — d v A^\ The equations of motion are 



N K N 



ein A® = 2 e fc „ Im W* n d^ n ] for k = 1, K . (A3) 
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Denoting by v n the VEV of ip n , the K x K mass matrix of the gauge fields is 

N 



Mm =2^u£e fc „e; n for k, I = 1, K . (A4) 



71 = 1 



The matrix M. is real and symmetric, so the system can be diagonalized with an orthogonal matrix U : U T M. U is 
diagonal with U^ 1 — U T . For the system to admit a massless gauge field, at least one eigenvalue of the mass matrix 
should vanish. Suppose that the j th eigenvalue \j vanishes, so that 

N 

X, = {U T M U).. =2j2 v n 4* = ( A5 ) 
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where we have defined 
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In the new basis 
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Ji|f ' is a massless candidate. It satisfies the equation 
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where F^J is its gauge field strength. There are two cases to consider. If all the scalar fields have a non-zero VEV, 



v n 7^ Vn = 1,...,N, then the condition (A5| implies that e 



that the massless gauge field A^P decouples from all the scalars, d v F$ = 



Vn 

Ci) 



1. 



., A. In this case, we see from (A8) 

On the other hand, for this field to 
remain coupled we need ej n ^ for at least one value of n. The condition ( A5 1 then implies that ip n has a zero VEV. 



However, since ij n ^ 0, we see from (IA8j) that Ap acquires an effective mass proportional to e jrl (\Sip n \ 2 ) due to its 
interaction with the fluctuations of ip n . Thus the U(l) gauge fields are either effectively massive or decoupled from 
the scalars and the other U(l) gauge fields. 



Appendix B: Lattice Formulation with 0(dx 2 ) Accuracy 

In this appendix we describe the discretized equations of motion that we evolve on the lattice. As is well known, 
special care has to be maid when discretizing a gauge theory in order to preserve gauge invariance on the lattice. 
For instance, this is necessary for constraint equations to follow from the dynamical equations that are evolved. The 
basic formalism that we use is standard (see e.g. |57| ) but, in the presence of gravity waves, we have to generalize it 
in order to reproduce the continuum theory up to 0(dx 2 ) and 0(dt 2 ) accuracy in the lattice spacing dx and timestep 
dt. 

In the lattice formulation, the scalar fields x an( i <P are defined at spacetime lattice points x while the gauge field 
A^ is defined in the middle of the segments between lattice points, at x + fi/2 where fi is a vector of length dx 11 in 
the /i direction, dx 1 = dx is the lattice spacing and dx° = dt is the timestep. In this appendix, we denote the lattice 
fields with a tilde to distinguish them from their continuum analogs. One then introduces link variables, also defined 
in the segments between lattice points, and related to the gauge field according to 



(Bl) 



In this equation the repeated index fj, is not summed. In the case of abelian symmetry, we can choose to treat either 
the gauge field or the links as the fundamental objects that we evolve numerically. We found that the first option 
led to faster simulations, so in the following we will express everything in terms of the gauge field. The discussion is 
easily generalized to the case where the links are considered as the fundamental variables. 

The starting point is the lattice action ( 40 1 which is a discretized version of the continuum action ([2| . The lattice 
expressions for the forward partial derivative, the forward gauge covariant derivative and the gauge field strength are 
given respectively by 



1 

dxt* 



d£x = ^ [x(x + A) - xO)] 



U a [x + 



(p{x + fi) — <p(x) 



Ffiv — ®u,A-v ~ &i Ap 



(B2) 
(B3) 
(B4) 



Again, repeated indices are not summed in Eqs.(B2 B3). In the limit dx^ — > 0, these expressions reduce to the 



continuum partial derivative, gauge covariant derivative and gauge field strength, respectively. The definitions (B2 



B4) imply that the action (40) is invariant under the lattice gauge transformation (41), which is a discretized version 



of a continuum gauge transformation. 
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The equations of motion following from ( 40 1 are obtained by the lattice equivalent of functional differentiation. 
They read 



d^F^ = 

where we have defined the backward partial derivative 

1 



dV_ 

dx 
dV 

dip* 

-2elm [if* D+p] 



" A ~ dx* ^ ~ ^ X ~ ^ 



and the backward gauge covariant derivative 



1 

dx* 



<p(x) — Ufilx 



<p(x- A) 



(B5) 

(B6) 
(B7) 

(B8) 
(B9) 



where again repeated indices are not summed. 

We evolve these equations in the temporal gauge (A = 0, Uq = 1) with the staggered leapfrog method, where the 
fields and their time derivatives are evaluated at times that differ by half a timestep. Explicitely, we have 



dt 



X\x + 



dt 

¥> [ x + y 



- / i dt 



= X \x 



dt 
dt 



dt 



~ I i dt 



di dfx - 
D^D+tp 



dV 
dx 
dV 

dip* 



(x) 



(x) 



together with f(x + dt) = f(x) + dt f (x + § ) and Ai 



Finally, the discretized version of Gauss constraint (|9j reads 

drEi = 2elm [p>* $\ 



dt dj Fji + 2e Im [ip* D+ ip] 
+ dt) = Ai (x + |) + dtEi (x + | + f V 



(BIO) 
(Bfl) 

(B12) 



(B13) 



ft follows from the dynamical equations (Bit Bf2) if it is satisfied initially, because these equations are derived from 
the same, gauge-invariant lattice action (40). Gauss constraint was satisfied down to machine precision in all our 
runs. 

In the limit dt — > 0, one can show that the equations of motion ( Bf OflBl 2 1 implies that the total energy density 



P : 



1 r 2 

2 X 



I ^ ^ 



1 z. 



- Ei E l + - 8+x d+x + DfCp (D+p) + - F %] F %] 



V 



(B14) 



is conserved, dp/dt = 0. Here (...) denotes the average over all the lattice points. Energy was conserved up to ~ 0.1% 
in all our runs. 

It is important to note that these discretized equations of motion reproduce the continuum ones up to 0(dt 2 ) and 
0(dx 2 ) when dt, dx —> 0. This can be checked explicitely, remembering that the gauge field and thus its equations of 
motion are defined in the segments between lattice points. In particular, we have 



DTDJCp ~ DiDnp(x) + 0{dx 2 ) 
£. u ( , i , 3 



*3 



d- Fjl 



*3 



djFji x - 



+ 0(dx 2 ) 
0(dx 2 ) 



Im [p* D+<p] ~ Im [p* A^ (x+ | } + 0(dx 2 ) 



(B15) 
(B16) 

(B17) 
(B18) 
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We now come to the discretized equations of motion for gravity waves hij(x), that we define at the lattice points 
as the scalar fields. Fortunately, we don't have to start from a discretized action for them because GW are gauge 
invariant and no constraint equation is associated to them. However, in the presence of gauge fields, special care is to 
be paid in order to reprod uce the continuum equation (31) up to 0(dx 2 ) and 0(dt 2 ) accuracy. In order to do so, we 



cannot replace the terms (|20|) in n^ T by their lattice ana 



ogs defined above, because these reduce to the continuum 
values at lattice points upTo O(dx) or O(dt) accuracy only. We thus have to construct new lattice expressions that 
reproduce the continuum up to second order and which lead to a gauge-invariant stress-energy tensor. For the partial 
derivative of the inflaton, this is just the symmetric derivative 



dfx 



1 

2dx 



X(x + i) - x(x - i) 



(B19) 



which indeed reduces to di\{x) 
the Higgs, by defining 



0(dx 2 ) for dx — > 0. We can achieve the same for the gauge covariant derivative of 




(p(x +i) — Ui 
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<p(x — i) 







(B20) 



Note also that, under the lattice gauge transformation (41 1, this transforms as Dfip — > e la Dfip, as it should 



For the lattice gauge field strength F«, note that it reduces with 0{dx ) accuracy to the continuum Fij evaluated 



in between lattice points, see Eq.(B16). To achieve the same at the lattice points themselves, we consider the clover 
average 




(B21) 



Similarly, the lattice electric field Ei reduces with 0(dx 2 ) and 0(dt 2 ) accuracy to its continuum conterpart evaluated 
at x + i/2 + dt/2. We thus consider its clover average 



E?(x) 



E, x 




EAx 




E, x 




Ei x 



(B22) 



which reduces to Ei(x) + 0{dx 2 ) + 0{dt 2 ) for dx, dt — > 0. In terms of these, the lattice analog of (20) is given by 

TT 



ri TT 



d?xdfx + 2Re 



(B23) 



It fulfills the conditions stated above, namely it is invariant under the lattice gauge transformation (41 ) and reproduces 
the continuum with second order accuracy in dx and dt. 

Contrary to Ei which was displaced by h alf-tim esteps, the clover average of the electric field Ef is defined at the 
timesteps themselves, as the other terms in (B23). This allows to preserve the leapfrog scheme for the gravity waves 



hij I X 



dt 



dt 



dt 



TT 



(x) 



(B24) 



This is easily implemented at each iteration, by first advancing Ei by half a time step, then advancing the gravity 
waves by a full timestep and finall y advancing Ei by the remaining half timestep. 

As discussed in subsection IV B the 0{dx 2 ) calculation of GW allows for a much better control on the UV part of 
the GW spectra. 
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